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PROJECT SUMMARY 

This project is motivated by the desire to develop numerical methods that will be useful 
in the study of compressible flows that exhibit aeroacoustic phenomena. Solutions to linear 
problems have been investigated through the development of a computer code based on the recent 
dispersion-relation-preserving (DRP) methodology. In regard to nonlinear problems, the class 
of essentially nonoscillatory (ENO) schemes have been considered as the primary candidates for 
solving aeroacoustic problems in which discontinuities are involved. Discontinuities in the solution 
itself ( e.g . shocks) as well as in the geometry on which the problem is defined have been studied. 
Two-dimensional nonlinear problems were considered in order to determine if the one-dimensional 
results obtained in the first phase of this project were extendable to more realistic problems. 
Conclusions have been drawn in regard to the ability to numerically predict solutions of nonlinear 
problems with shocks to high-order accuracy. 
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1. Report Summary 

In the remainder of this report, an overview of each of the technical problems and a synopsis 
of their achievements are presented in Section 2 and Section 3. The details are available in the 
technical publications listed in Section 4 and are sampled in Section 5 

2. High-Order Accurate Methods for Compressible Flows with Shocks 

Phase 1 Research 

During the first phase of this project, reported in (Casper and Baysal, 1995), a prime 
objective was to determine the requirements of a numerical to predict a high-order accurate 
solution of a compressible with shocks. This subsection provides an overview of this Phase 1 
work. 

A fourth-order accurate ENO-based code was used to investigate the solutions of inviscid, 
compressible flows with shocks. Design accuracy was achieved in the smooth regions of a steady- 
state, quasi-one-dimensional, shocked nozzle flow. However, in an unsteady one-dimensional test 
case, only first-order results were obtained downstream of a sound-shock interaction. Similar 
results were obtained with a fourth-order compact finite-difference code. The difficulty in 
obtaining a globally high-order accurate solution to an unsteady discontinuous problem with a 
shock-capturing method was examined through the study of a simplified, linear model problem. 
The ENO code was modified with a subcell resolution technique in order to obtain a globally 
high-order accurate solution. The sound-shock interaction problem was again solved with the 
modified scheme, and a globally fourth-order accurate solution was achieved. Details of this 
research can be found in (Casper and Carpenter, 1995). 

In addition to questions concerning the accuracy of discontinuous flows, the issue of 
discontinuous geometries was also a topic of interest. The performance of a fourth-order 
ENO method was tested on smooth solutions defined on piecewise smooth geometries. The 
propagation of sound in a quasi-one-dimensional nozzle was considered as a test case. Not only 
was the geometry for this problem not smooth to the order of the numerical scheme, the mesh 
was a discontinuous function of the nozzle length. The code was modified through the adaptive 
stencil strategy employed by ENO methods. The solution’s accuracy, with respect to entropy, 
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was found to be fourth-order accurate, as per the design of the numerical scheme. Details of this 
research can be found in (Casper, 1995). 

Phase 2 Research 

The difficult issues raised by the results obtained in Phase 1 are topics of the second phase 
of this project and are discussed below. The sound-shock interaction was revisited, in order to 
study the relative merits of linear and nonlinear methods, as well as to quantify the behavior 
exhibited by these methods in the presence of solution discontinuities. The issue of extension of 
these methods to more than one spatial dimension was also of concern. Because the extension 
of subcell resolution techniques to multiple spatial dimensions is impractical, the relative merits 
of high-order shock-capturing methods had to be reconsidered. A two-dimensional problem 
was solved in order to compare various methods for their accuracy with respect to their shock 
resolution. 

A sixth-order compact finite-difference scheme was used to investigate various one-dimensional, 
discontinuous flows: steady-state flow in a nozzle, sound-shock interaction, time-dependent 
Burgers’ equation, and a linear model problem. The results obtained from the solution of the 
Euler equations reaffirmed the conclusions from the previous research reported in (Casper and 
Carpenter, 1995). It was found that, although both the linear compact scheme and the fourth- 
order ENO method were high-order accurate in the steady Euler case, the ENO results were 
significantly better, perhaps owing to the ability to capture shock in a nonoscillatory fashion. 

For both schemes, the unsteady Euler solutions were only first-order accurate downstream of 
the shock. However, the solution of the unsteady Burgers’ equation was found to be of design 
accuracy away from the discontinuity. An interesting result was obtained by revisiting the linear 
model problem with discontinuous wave speed. It was found that the first-order error made across 
the discontinuity was largely a phase error, and that post-processing would enable second-order 
accuracy to be recovered. Details of this research can be found in (Carpenter and Casper, 1996). 

A fourth-order compact finite-difference scheme and a third-order ENO scheme were used 
to predict the two-dimensional, steady flow around a supersonic cylindrical projectile. Because 
of the previous research that obtained high-order results for a steady one-dimensional shocked 
flow, high-order results might have been expected for the two-dimensional steady case. However, 
high-order accuracy was obtained with the ENO method only when the shock was aligned with the 
computational mesh, in which case the problem becomes one-dimensional. When the problem 
is determined such that the shock is skew to the mesh, only first-order accuracy is obtained. It 
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was concluded that such errors are inherent in numerical shock-capturing methods, as presently 
understood. A better understanding of high-order accurate shock-capturing methods in a general 
multi-dimensional context will be required before CFD can play a significant role in aeroacoustic 
applications that involve flows with shocks. Details of this research can be found in (Carpenter 
and Casper, 1997). 

Copies of all reports pertaining to the results obtained from the research associated with 
this project are included in Section 5. 

3. Low-Dispersion Schemes for the Direct Computation of Acoustic Wave Propagati on 

Computational aeroacoustics (CAA) may be defined as the application of numerical tech- 
niques for the direct calculation of aerodynamic sound generation and propagation starting from 
the first principles. Most computational fluid dynamics (CFD) schemes, however, are not ad- 
equately accurate for solving the aeroacoustics problems (Lighthill, 1992). Their amplitudes 
are often orders of magnitude smaller, and yet the frequencies are orders of magnitude larger 
than the flow field variations generating the sound. Further, high-fidelity is paramount for the 
resolution of acoustic problems; but a consistent, stable, and convergent, high-order scheme is 
not necessarily dispersion-relation preserving and thus does not necessarily guarantee a good 
quality numerical wave solution for an acoustic problem. Hence, among the requirements that 
should be placed on a CAA algorithm are the minim al dispersion and dissipation features. 

Therefore, a fourth-order accurate dispersion-relation-preserving (DRP) method (Tam and 
Webb, 1993) was previously developed for a variety of wave propagation problems by solving 
the linearized Euler equations. Evaluation of the developed baseline algorithm and the spectral 
analysis methods for time-series data were reported in a paper (Vanel and Baysal, 1995). In this 
investigation, a number of observations and recommendations were made for future investigations 
to extend the scheme to solve some application problems. The present investigation started 
precisely with this impetus. 

The major point of the work accomplished in the first year was a robust low-dispersion- 
method for the linear Euler equations which was also relatively more versatile and efficient: 

(i) The previously written computer code (Baysal, 1995a, Vanel and Baysal, 1995) for the 2D 
DRP algorithm has been completely revised for a number of reasons. The new version is written 
using dynamic memory allocations, for example, and it is more tailored to run on workstations. 


4 



(ii) The boundary condition routines for the linear DRP method (radiation and outflow boundary 
conditions) have been complemented for solid walls. 

(iii) For improved wave propagation simulations, selective artificial damping has been incorpo- 
rated and thoroughly tested. The user can now select the band of the error that is wished to be 
suppressed. 

(iv) The method has been tested for various types initial-value and periodic problems generating 
incident and reflected waves off the walls either inside of the domain or on a side of the domain 
(from The First Workshop on Benchmark Problems, Hardin et al., 1995). 

(v) The method has been extended and the optimized coefficients have been derived for a spatially 
sixth-order accurate scheme. The spurious waves and nonisotropicity has been reduced to about 
machine zero. 

(vi) The linearized diffusion terms have been derived and discretized with the DRP scheme and 
the scheme has been extended for the linearized 2D Navier-Stokes equations. 

(vii) . An optimized, variable-order Runge-Kutta method has been implemented as a low-storage 
alternative to the time integration method. 

(viii) The extension of the method has been successful for the generalized curvilinear coordinates. 

Several issues were investigated. Among them were higher-order accuracy, choice of 
boundary conditions and differencing stencils, effects of viscosity, low-storage time integration, 
generalized curvilinear coordinates, periodic sources, their reflections and interference patterns 
from a flat wall and scattering from a circular cylinder. The results were found to be promising 
en route to the aeroacoustic simulations of realistic engineering problems. 

In the second year of the project, the major points of achievement were solving the scattering 
problems given in the Second Workshop on Benchmark Problems (Hardin and Tam, 1996, 
1997), and the nonlinear wave propagation. By and large, the simulations of the propagation 
of acoustic waves, their reflections and scattering, in particular, the initial-value problem with 
the acoustic pulse, were successful (Bay sal and Kaushik, 1997). Two necessary building blocks 
to success, once a suitable CAA scheme was selected, were the correct boundary conditions, 
and an adequate arid efficient grid. Notwithstanding the imperfectly orthogonal grids and 
the required transformation metrics, employing the body-fitted coordinates allowed a straight 
forward implementation of the boundary conditions. 
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There axe certain important applications, e.g. the jet screech in a wake-vortex interaction 
flow, boundary layer noise, cavity noise, etc., where the assumption of linear wave propagation in 
uniform flow would be too compromising. Therefore, the linear dispersion-relation-preserving 
scheme and its boundary conditions were extended for the nonlinear Euler equations (Baysal 
et al., 1997). This allows computing, a nonuniform flowfield and a nonlinear acoustic wave 
propagation in such a medium, by the same scheme. By casting all the equations, boundary 
conditions, and the solution scheme in generalized curvilinear coordinates, the solutions were 
made possible for non-Cartesian domains and, for the better deployment of the grid points, 
nonuniform grid step sizes are also possible. 

The methodology has been tested for a number of simple initial- value and periodic-source 
problems. The wall boundary condition derived from the momentum equations and implemented 
through a pressure at a ghost point, and the radiation boundary condition derived from the 
asymptotic solution to the Euler equations, were shown to be effective for the nonlinear equations 
and nonuniform flows. The non-reflective characteristic boundary conditions had success limited 
to the nonlinear waves in no mean flow, and failed for nonlinear waves in nonuniform flow. 

The methods are now being used to simulate the following two "simplified application" 
problems: a) noise generated by a 2D, steady, supersonic jet flow, b) noise generated by a 2D 
cavity flow. 

Finally, copies of the papers reporting the results of this project are enclosed in Section 5. 


Principal Investigator: Dr. Jay Casper 
Co-Principal Investigator: Dr. Oktay Baysal 
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LOW-DISPERSION SCHEME FOR NONLINEAR ACOUSTIC WAVES 
IN NONUNIFORM MEAN FLOW 


Oktay Baysal,* Dinesh K. Kaushik,** and Moumen Idres** 


Aerospace Engineering Department 
Old Dominion University 
Norfolk , Virginia 23529-0247 USA 

The linear dispersion- relation-preserving scheme and its boundary conditions have been extended to the 
nonlinear Euler equations . This allowed computing , a nonuniform flowfield and a nonlinear acoustic wave 
propagation in such a medium , by the same scheme. By casting all the equations, boundary conditions , and the 
solution scheme in generalized curvilinear coordinates, the solutions were made possible for non-Cartesian 
domains and , for the better deployment of the grid points , nonuniform grid step sizes could be used. It has 
been tested for a number of simple initial-value and periodic-source problems. A simple demonstration of the 
difference between a linear and nonlinear propagation was conducted. The wall boundary condition , derived 
from the momentum equations and implemented through a pressure at a ghost point , and the radiation boundary 
condition , derived from the asymptotic solution to the Euler equations , have proven to be effective for the 
nonlinear equations and nonuniform flows. The nonreflective characteristic boundary conditions also have 
shown success but limited to the nonlinear waves in no mean flow , and failed for nonlinear waves in 
nonuniform flow. 


Introduction 1 

There are certain aeroacoustics applications, e.g. 
the jet screech in a wake- vortex interaction flow, 
boundary layer noise, cavity noise, where the 
assumption of linear wave propagation in uniform flow 
would be too compromising. Consequently, there were 
some recent progress in the attempts to compute the 
propagation of nonlinear acoustic waves in a 
nonuniform flowfield using high-order schemes. For 
example, Tam and Shen 1 have reported some extensions 
of a dispersion-relation-preserving scheme 2 in this 
direction and solved the inviscid Burger's equation. 

Considering accuracy as well as efficiency, some 
success has so far been reported in direct computations, 
but they are predominantly limited to modeled periodic 
noise generation and propagation by assuming linear 
waves in uniform flow.^’ 4 In addition to the cases 
reported in the First and Second CAA Workshops on 
Benchmark Problems,^ 4 the authors have shown some 
successful applications with low-dispersion schemes by 
solving the linearized Euler or Navier-Stokes 
equations.^" 7 However, these equations would not be 
able to represent an aeroacoustic phenomenon, such as, 
that observed in the case of a high-speed flow past a 
cavity** or the engine internal noised One way to 
approach such problems would require the nonlinear 
equations, and the extension of low-dispersion, high- 
order schemes and boundary conditions for these 
equations. This is precisely the goal of the present 
paper, albeit a modest step in that direction. 


i* 

^Professor and Eminent Scholar. Associate Fellow, AIAA. 
Graduate Research Assistant. 

Copyright © 1997 by the American Institute of Aeronautics and 
Astronautics, Inc. All rights reserved. 


Methodology 

The present model employs the two-dimensional, 
compressible, nonlinear Euler equations. Their 
nondimensionalized and nonconservation form is derived 
in the generalized curvilinear coordinates. It is 
contended that, unless strong gradients are expected in 
the nonuniform mean flow, there should not be 
discernible differences in employing the equations in 
conservation or non-conservation forms 1 : 
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The density (p), pressure (p), and velocity (w,v) are 
normalized using the freestream values of speed of sound 
( ) and density. The length is normalized either by 
the step size (Ax) or the cylinder diameter (D). The 
metric terms ( % x , etc.) are computed analytically for the 
simple grid cases considered herein. 
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For the spatial integration, the dispersion-relation- 
preserving scheme, 2 which was demonstrated for various 
linear cases, 5 " 7 is now extended for the nonlinear 
equations in generalized coordinates. Even with second- 
order schemes, 8 when used in resolving the highly 
nonlinear flow or acoustic phenomena, it is common to 
have limiters or artificial damping mechanisms, either 
implicit in the scheme (upwind schemes) or explicitly 
added (central schemes). Since simulating low- 
amplitude acoustic waves for a long time and for many 
wavelengths of travel is the objective in CAA, devising 
such artificial damping mechanisms requires extreme 
care. Constant damping over all wave numbers need to 
be avoided; rather, selective damping, 10 which has 
small damping over the long wave range but significant 
damping in the short wave range, should be devised. 

The present central-difference scheme includes similar 
terms ( D^ m ) to overcome the expected spurious 
oscillations in a highly nonlinear region: 

(#)".» = -Km + Si* + D lm > ( 4 ) 

where 

M , M 

Rlm-felElm JpQf +J , m ]+^[Fl m JpQl m+j ] (5) 
and 

D lm .Zdjl-^p-Qe+j.m +^T Qt.m+jb Re D =' 2 ^' 

( 6 ) 

In eqs. (4) to (6), l and m are the spatial indices and n 
indicates the time level. The coefficients aj and dj were 
computed by minimizing the dispersion-relation error. 5 " 
7 For N=M , difference is central, for N=0, it is fully 
forward, and for M- 0, it is fully backward. However, 
since these high-order stencils require multiple layers of 
boundary cells, all combinations between a central and a 
fully-one-sided difference need also be derived. Only 
then, it would be possible to always utilize the 
information from the nearest possible points for better 
accuracy. In the present computations, a fourth-order 
scheme was used, requiring 7-point stencils: N takes 
values from 6 to 0, M takes values from 0 to 6, and 
N+M=6. The interior ceils were computed using central 
differences with JV=M=3. 

For the temporal integration, a low-storage, low- 
dispersion JRunge-Kutta scheme, 11 which was 
previously,^ 7 demonstrated for various linear cases, is 
now extended for the nonlinear equations in generalized 
coordinates: 

Q (°) m Qn t Q U) = q(0 ) 4f( ij2)(«W) > 

Q n+l =Q W , i = l,2,..,p (7) 


Ci=UPp-k+2’ i = 2,...,p . (8) 

k—2 

The coefficients c[ were computed by considering the 
amplification factor of the Runge-Kutta scheme, then 
minimizing the dispersion-relation error. The time 
steps were determined from the stability as well as the 
accuracy limits. In the present study, a five-stage (p =5) 
Runge-Kutta was used, which required two levels of 
storage and it was at least second-order accurate. 

One of the challenges en route to the objective was 
clearly the implementation of the boundary conditions. 
Although a number of papers have already reported 
various types of nonreflecting or absorbing boundary 
conditions, in representing all forms of waves in 
nonuniform mean flow, even at the boundary, few were 
derived for the Euler equations. The boundary 
conditions proposed by Tam and Dong 12 also assume 
weakly nonuniform flow at the boundary, hence can be 
based on the asymptotic solutions. Further, it requires 
some knowledge of the total mass flux in constructing 
the mean flow quantities. An improvement by Dong 
and Mankbadi 9 has eliminated the latter but not the 
former impediment. Further analysis of this proposed 
improvement is reported by Dong. 13 The present 
method implemented these boundary conditions. 

Radiation boundary conditions are derived from the 
asymptotic solutions, 12 ’ 13 assuming the mean flow is 
only weakly nonhomogeneous and, in the farfield, the 
outgoing disturbances are propagating in the radial 
direction relative to the noise source at speed V w . 
Denoting the acoustic field properties as the difference 
between the total ( Q ) and the mean ( Q ) properties, they 
are: 


^ +A m^ +B *^21+C(Q-Q) = 0 (9) 


where 


A = Vw f^,B = V w ^^, C = ^,r = ^7+y 2 


( 10 ) 


V w =u cos0 + vsin0 + -ya 2 -(« sin£~ vcos0) 2 (1 1) 


In addition to the acoustic waves, however, 
vorticity and entropy waves also reach the downstream 
boundary. For the pressure, the outflow boundary 
condition, 12 is still derived from eq. (9). However, the 
density and velocity, lumped into a vector, Q\ are 
derived to be as follows: 


^+M(Q'-~Q') = rhs , where (12) 


and 


^ * v i V 1 >- 

-W,*5 a +i, i! £ a n' 


(13) 


The indices p and i indicate the order and the stage of the 
method, respectively. As for the coefficients, pj =0 and 
the other coefficients pi were determined from, 


2 



Low-dispersion scheme for nonlinear acoustic waves 
in nonuniform mean flow. 

Baysal, Kaushik, Jdres 


AIAA 97-1 582 CP 


To maintain the impermeability at , for example, 
an T) = 0 solid wall, the normal contravariant velocity 
is set to zero, V = rj x u + rj y v = 0 . Using this condition, 

an equation for the normal pressure gradient is obtained 
from the two momentum equations, 




By DRP-type finite differencing, 14 this equation at time 
level n and wall location (7,>v), provides the ghost 
point value (subscript -1) of pressure: 




Pe -- 1= t £o ajPtj] ~ 

(pU), , \ M / v M 

. + , /1C ^ 

+ 25 Z,OjPt+j tW } . ( 15 ) 


In order to derive the nonreflecting characteristic 
boundary conditions as an alternative to all of the above, 
first a locally one-dimensional ( l ,i = 7,2) flow is 
assumed, then the Euler equations are written in their 
characteristic form, 15 


^ T P a ^ +r ‘u z:0 ' ^r~° 2 i^ +r 2~ 0 ( 16 ) 

where 

r 2 - v ‘^-° 2 f 1 (n) 

i = 7,2 {no sum) 


In order to have a nonreflecting condition, 15 ’ 16 at a 
given grid point on a boundary, the local perturbations 
propagated along incoming (to computational domain) 
characteristics are made to vanish, and perturbations 
along outgoing characteristics are computed by the one- 
sided DRP differences. This requires checking the local 

Mach number ( M l s ■&-), and the signs of the directed 


unit normal ( n l 



no sum) and its scalar product 


with the normal contravariant velocity 
( U l = X ihh ). The tests for the proper choice at a 

given boundary point in the computational plane are 
summarized in table 1. 


Results 

Previously, 5 ' 7 the linear algorithm was developed 
and demonstrated for the following test cases in uniform 
flow: (a) single acoustic pulse, (b) three successive 
acoustic pulses, (c) two simultaneous acoustic pulses, 
(d) acoustic pulse near a flat wall, (e) periodic source 
near a flat wall with and without a circular bump, (f) 


acoustic pulse near a circular cylinder (scattering 
problem), (g) periodic source near a circular cylinder 
(scattering problem). The present nonlinear algorithm 
and its boundary conditions are now being evaluated by 
considering the following numerical examples (table 2): 
(1) acoustic pulse generated wave propagation computed 
by the nonlinear algorithm contrasted with that by the 
linear algorithm; 5 ' 7 (2) acoustic pulse reflected from a 
flat plate computed with the nonreflective characteristic 
boundary conditions compared with the solution using 
radiation boundaiy conditions; (3) acoustic pulse 
reflected from a bump on a flat plate; and (4) periodic 
acoustic source in a potential sink ( nonuniform ) flow. 

The acoustic waves were generated either by a 
periodic source to generate both the nonuniform mean 
flow and the acoustic waves, 

s(x,y,t) = s + s' coscot (18) 

where 

J = -eexp{-a[U~^) 2 +(y-~y iS ) 2 ]} (19a) 

s' -e exp {- <7 [(x -x s ) 2 +(y- y s ) 2 ]} (1 9b) 

am^ 9 am Jgf . (19c) 

or by an initial acoustic pulse introduced to the field at t 
=0 by setting u =v =0 and s(x,y) = s ' . In eq. (18), the 
frequency is co =0 . 2rc and the rest of the parameter 
values are given in table 2 for each case. 

The propagation in a no mean flow medium of the 
acoustic waves generated by the pulse (case 1), are 
presented in fig. 1. When computed by the linear 
equations, 5 ' 7 the waves travel at a constant speed, hence 
preserving their shape. Their convection velocity, when 
computed by the nonlinear equations (1), is a function 
of local state properties, which in turn, are functions of 
space and time. This results in a successor wave point 
taking over the predecessor wave point, manifesting in 
an ever steepening wave front (fig. 1). For case 2, a 
solid wall is added as the lower boundaiy of the domain 
(figs. 2 and 3). The interference pattern of the incident 
and the reflected waves is computed by the nonlinear 
algorithm. Until the waves reach the non-solid 
boundaries, there virtually is no difference between the 
two solutions. From t =54 and on, minor differences 
between the nonreflective characteristic boundary 
conditions (fig. 2 and eq. (16)) and the radiation 
boundary conditions (fig. 3 and eq. (9)), are visible. 

To generate reflected waves that are nonsimilar to 
the incident waves, and to check the success of the 
algorithm with non-rectangular boundaries, hence 
requiring non-Cartesian treatments, such as, the 
curvilinear coordinates employed herein, in case 3, a 
circular bump is placed on the flat plate (fig. 4). The 
nonlinear wave propagation and the interference are 
successfully computed with the present algorithm and 
the radiation boundary condition. 

In case 4, a nonuniform flow is generated by 
placing a potential sink (s in eq. (18)) and computed by 
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the present nonlinear algorithm. The residual is driven 
to the machine zero (fig. 5a), and the computed 
distributions of density and velocity along the x-axis are 
successfully compared with the exact solution 12 (fig. 

5). Then, the acoustic waves generated by the periodic 
source are computed and their nonlinear propagation 
(distributions of p acoustic =P~P and V acoustic = V- V ) 
are compared (fig. 6) with the exact solution. 12 ' 13 

Conclusions 

The linear dispersion-relation-preserving scheme 
and its boundary conditions have been extended to the 
nonlinear Euler equations. This allowed computing, a 
nonuniform flowfield and a nonlinear acoustic wave 
propagation in such a medium, by the same scheme. 

By casting ail the equations, boundary conditions, and 
the solution scheme in generalized curvilinear 
coordinates, the solutions were made possible for non- 
Cartesian domains and, for the better deployment of the 
grid points, nonuniform grid step sizes could be used. 

The method has been tested for a number of simple 
initial-value and periodic-source problems. A simple 
demonstration of the difference between a linear and 
nonlinear propagation was conducted. The wall 
boundary condition derived from the momentum 
equations and implemented through a pressure at a ghost 
point, and the radiation boundary condition derived from 
the asymptotic solution to the Euler equations, have 
proven to be effective for the nonlinear equations and 
nonuniform flows. The nonreflective characteristic 
boundary conditions also have shown success but 
limited to the nonlinear waves in no mean flow, and 
failed for nonlinear waves in nonuniform flow. 

Currently, work is in progress to further improve 
the boundary conditions. Also, the scheme is being 
extended to improve its computational efficiency. 
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Tables and Figures 


Tabl e 1 . Nonreflective characteristic boundary conditions. 
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£/' •«' <0 

M‘ <1 
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O 

r;,r 2 by 

forward DRP 

Jj set to zero 

n by 

forward DRP 

r 2 >r 3 

set to zero 

Si 

V 

O 

T2.-T? by 
backward DRP 

17 set to zero 

O by 

backward DRP 

n,r 2 

set to zero 

M‘ > 1 

A 

O 

P 1 ^ 2 X 3 by forward DRP 

rj,r 2 ,r 3 set to zero 

V 

O 

ri,r 2 ,r 3 by backward DRP 


Table 2. Description of computational cases. 


Case 

£,£ 

b y b 


Grid type, 
size 

At, At 

Re 

Figure 


0, 1 

0, 6Ax 

0,0 

H 141x71 

0, 3.7E-1 

1E1 

1 

2: pulse reflected from 
flat plate; no mean flow 

0,1 

0, 6Ax 

0, 20 Ax 

H 141x71 

0, 3.7E-1 

1E1 

2,3 

3: pulse reflected from 
bump on flat plate; 
no mean flow 

o, 1 

0, 0.6D 

0, 2D 

O 65x65 

0, 8E-3 

1E3 

■ 

4: periodic source in 
nonuniform potential 
sink flow 

6.128E-2, 

IE-3 

6Ax, 

3Ax 

0,0 


2E-1, 

8E-2 

1E1 

5,6 



Fig. 1 Linear and nonlinear wave propagation (case 1). 
Time history of the positive part of the wave. 
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Fig. 2 Reflection of a 
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t=4.8 t=5.6 


Fig. 4 Reflection of an acoustic pulse from a bump on a flat plate (case 3). Acoustic 
pressure contours computed using asymptotic boundary conditions 


log10(max(dQ/dt)) 



(a) 




Fig. 5 Nonuniform mean sink flow (case 4) . (a) Residual history, (b) Density distribution 
and (c) Velocity distribution , along x-axis. 




Fig. 6 Nonlinear acoustic wave propagation from a periodic source (case 4). Instantaneous 
(a) Density distribution, (b) Velocity distribution, along x-axis. 
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COMPUTATION OF ACOUSTIC SCATTERING BY A LOW-DISPERSION 

SCHEME 


Oktay Baysal and Dinesh K. Kaushik 

Department of Aerospace Engineering 
Old Dominion University, Norfolk, Virginia 23529-0247 


Physical problem and background 

The objective is the evaluation of a proposed computational aeroacoustics (CAA) method in simulating an 
acoustic scattering problem. An example may be the sound field generated by a propeller scattered off by 
the fuselage of an aircraft. The pressure loading on the fuselage would be an input to the interior noise 
problem. To idealize the problem, the fuselage is assumed to be a circular cylinder and the noise 
generation by the propeller is represented by a line source. 

A typical CFD algorithm may not be adequate for this aeroacoustics problem: amplitudes are an order 
of magnitude i? smaller yet frequencies are d larger than the flow variations generating the sound. For 
instance, an 2) CFD method was previously used for a nonlinear wave propagation problem in 
unsteady, nonuniform mean flow (Baysal et al., 1994). It was observed that a direct simulation of acoustic 
waves using a higher-order CFD would become prohibitively expensive, due to the required excessive 
number of grid points per wavelength (PPM). Also, CAA would need minimal dispersion and dissipation, 
which would preclude a typical t? (2) CFD method from a long-term wave propagation simulation. 
Furthermore, a consistent, stable, convergent, high-order scheme is not necessarily dispersion-relation 
preserving, i.e. no guarantee for a quality numerical solution. Therefore, a baseline i?( 4) dispersion- 
relation-preserving (DRP) method (Tam and Webb, 1993) was investigated (Vanel and Baysal, 1997) for a 
variety of wave propagation problems, such as, single-, simultaneous-, and successive-acoustic -pulses. 
Then, a number of algorithmic extensions were performed (Kaushik and Baysal, 1996), when the 
following were studied: viscous effects by solving the linearized Navier-Stokes equations, low-storage 
and low-CPU time integration by an optimized Runge-Kutta scheme, generalized curvilinear coordinates 
for curved boundaries, higher-order accuracy by comparing d(4) DRP vs. d(6) DRP, and choice of 
boundary conditions and differencing stencils. The scheme is now being investigated for nonlinear wave 
propagation in nonuniform flow (Baysal et al., 1997). 

Mathematical approach 

The linearized, two-dimensional, compressible, Euler equations were considered in generalized curvilinear 
coordinates 

(1) <E = -R(U) + S, where R(U) = ^+^- . 

dt d£ drj 

The primitive variables, U, and the transformed fluxes, E and F, are, 

(2) U = i[p u v pf , E=i[Z x E+l; y F}, F=jr[ri x E + rj y F] . 
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and the physical fluxes and the source vector are. 


u + M 0 p 


V 


'O' 

M 0 u + p 


0 


0 

M 0 v 

, F = 

P 

, 5 = 

0 

MqP + u 


V 


.5 ' 4 . 


The perturbed values of density (p), pressure (p), and velocity (u,v) are denoted without a subscript, but 
those of the mean flow are demarcated using the subscript 0. These variables were normalized using the 
cylinder diameter (D) for length, speed of sound ( ao ) for velocity, D/ao for time, po for density, and 
Po( a o) 2 f° r pressure. 

Dispersion relation of a proposed numerical scheme should match closely that of eq. (1) for large 
range of resolution; i.e. a and a be close approximations to a and co. Assuming a Ax was a periodic 

(27t) function of aAx with, Fourier-Laplace transforms rendered. 


( 4 ) 


„ -i M 
a = — X a, e 
Ax j=-L J 


rajAx 


and 


® = - 4 ? rrr 

At I b i e v,mA ‘ 

j=0 


Finite difference coefficients were obtained from Taylor series expansions as one-parameter family, and 
the remaining coefficient from an error minimization, where the integration limit e depended on the 
shortest wavelength to be simulated: 


( 5 ) 


e 

min Ej = min J | aAx — a Ax\ d(aAx) 
-£ 


where 


\aAx\ -< e and e = 



for 


[4.5Ac cl 
{ 7 Ax J ' 


The time integration of eq. (1) was performed in two different ways. In the first approach, a four- 
point finite difference, which in a standard sense could be up to third-order accurate, was derived from the 
Taylor series as a one-parameter family. The remaining coefficient (bj) was determined, as in the spatial 
coefficients, by minimizing the discrepancy between the effective and the exact dispersion relations. After 

discretizing all the terms in eqs. (l)-(3), the resulting 0(At 2 ,Ax N+M ~ 2 ) DRP scheme was as follows: 


( 6 ) 


u: 


TJ^U^+AtlbjRti 

j—0 


M 

where R n t , m =^ I E n t+j m 

j=-N 


j=-N 


In eq. (6), i and m are the spatial indices and n indicates the time level. For N=M, difference is central, 
for N=0, it is fully forward, and for M=0, it is folly backward. All the interior cells were computed using 
central differences with N=M=3. However, since these high-order stencils require multiple layers of 
boundary cells, all combinations between a central and a folly-one-sided difference need also be derived. 
Only then, it would be possible to always utilize the information from the nearest possible points for better 
accuracy. In the present computations a fourth-order scheme was used, requiring 7-point stencil: N takes 
values from 6 to 0, M takes values from 0 to 6, and N+M=6. 

The numerically stable maximum time step At was calculated from the Courant-Friedrichs-Lewy 
relation. For example, for the fourth-order scheme in Cartesian coordinates, the stable CFL number was 
found to be 0.4. However, after analyzing the numerical damping of the time integration scheme, the CFL 
value was set to the more stringent value of 0.19. Since, however, time integration with DRP would 
require the storage of four time levels, a lower storage alternative, the low-dissipation and low-dispersion 
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five-stage Runge-Kutta scheme (Hu et al., 1996) was adapted and implemented. The resulting scheme had 
the spatial integration identical to eq. (6), but the time integration was replaced by the following: 


(7) U {0) = U n , U 0) =U (0) -/?,• AtR (i ~ n , U n+, = U {p \ i = 1,2,.., p . 

The indices n, p and i indicate the time level, and the order and the stage of the Runge-Kutta method, 
respectively. As for the coefficients, /3/=0 and the other coefficients A were determined from 

(8) c,- = FI Pp-k+2 ’ i = 2,...,p , 

k—2 

The coefficients c,- were computed by considering the amplification factor of the Runge-Kutta scheme, then 
minimizing the dispersion-relation error. The time steps were determined from the stability as well as the 
accuracy limits. In the present study, a five-stage Runge-Kutta (p =5) was used, which required two 
levels of storage and it was at least second-order accurate. When it was used with the present 7-point 
spatial stencil to solve the scalar linear wave equation, the CFL limit from the stability was found to be 
3.05, but it was only 1.16 from the accuracy limit. Since, however, this still was larger than the CFL limit 
of the DRP time integration, this method was also more efficient in processing time. 

Usually low-order schemes are used in resolving the highly nonlinear flow or acoustic phenomena. 
Even with second-order schemes, it is common to have limiters or artificial damping mechanisms, either 
implicit in the scheme or explicitly added. Since simulating low-amplitude acoustic waves for a long time 
and for many wavelengths of travel is the objective in CAA, devising such artificial damping mechanisms 
requires extreme care. For example, constant damping over all wave numbers need to be avoided. Tam et 
al. (1993) suggested terms, which have small damping over the long wave range but significant damping in 
the short wave range. The present central-difference scheme includes similar terms with a user specified 
artificial Reynolds number, to overcome the expected spurious oscillations: 


(9) 


L 

u l,m ~ R e . 


N 

I Cj[ 
j=-N 


AS 


2 U t+j,m + 


Ar] 


2 Ulm + j 1 , 


where 


_ Po a oD 


The boundaries should be transparent to the acoustic disturbances reaching them to avoid any 
degradation of the numerical solution. From the asymptotic solutions of the finite difference form of the 
Euler equations, a set of radiation boundary conditions, were derived and implemented. Therefore, 
following Tam and Webb (1993) and from the asymptotic solutions of the finite difference form of eq. (1), 
a set of radiation boundary conditions were derived, 


( 10 ) 


du . du D dU - n , 

+ A-—+B——+CU = 0, where 
dt dc, dr/ 


( 11 ) 


BmV 2^L t c =Jf and r =^x 2 +y 2 ; V = f M 0+A //-(A/ 0 ^ 


For an inviscid solid wall (w) or a symmetry plane, the impermeability condition dictates that the normal 
contravariant velocity be zero at all times: 


(12) v = ri x ii+ri y v = 0 . 

When this equation was constructed from the p- and ^-momentum equations, after multiplying them by the 
appropriate metrics and adding, the wall value of pressure at the ghost point (subscript -l) was obtained. 
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Pi-i = 


—^77 


Ad ^ ( \ ^ 

)^"Z f^ x ^ x h,w j^_p Ue+ j’ w + \ T] y^ x )i, w j]i^/ Vf+J ’ 


N+M {rht*+r, y Z y ) ( ' W 3 


N+M 


( 13 ) 


+ { r lx)i w Z a j u lJ ] + 'L a jPi+j,w} a [ % a jPt,j] 

j=w ’ ;'=-J ’ j=w 


The coefficients for all the boundary conditions (Tam and Dong, 1994) were derived by an analogous 
method to that of the boundary region cells. At the comers, two separate ghost points were used, one for 
each boundary, hence both boundaries’ conditions were satisfied. 

Results 

The present method and its boundary conditions were evaluated by considering a number of reflection or 
scattering cases (table 1), all in quiescent medium, i.e. Mq=0. The acoustic waves were generated either 
by an initial acoustic pulse introduced to the field at t — 0 by setting S =u =v =0, and 


(14) P- P exp {-ln2 [— — — ^ (Category 1_Problem 2), 

b 

or by a periodic source with co=8k. 


(15) S 4 = s exp{-ln2 [— — — ^-]} sin(fflt) (Category 1_Problem 1). 

b 


Table 1 . Description of computational cases 


Case 

Equation 

p or s 

b 

* s .ys 

Grid 

At 

Figure 

1 

14 

1.0 

3D 

4D, 0 

H 

251x101 

1 . 00 E -2 

1 

2 

14 

1.0 

0.2D 

4D, 0 

O 

401x181 

5 . 00 E -3 

- 

3* 

14 

1.0 

0.2D 

4D, 0 

O 

801x181 

5 . 00 E -3 

m mm 

4 

15 

0.01 

3D 

0, 2D 

H 

251x101 

1 . 00 E -2 

5 

5 

15 

1.0 

3D 

0, 2D 

H 

251x101 

1.00E-2 

6 

6 

15 

1.0 

0.2D 

4D, 0 

O 

801x181 

2.5E-3 

7 

Hi 

15 

1.0 

0.2D 

4D, 0 

O 

361x321 

2.50E-3 

1.25E-3 

8, 9 


Category l_Problem 2 Category l_Problem 1 


Gaussian pulse: an initial-value problem 

In case 1, a coarse H-grid was generated, where the £ -lines were along the cylinder and the centerline, 
and the rj -lines were perpendicular to the centerline. The transformed computational domain was 
rectangular with uniform steps in each direction and orthogonal grid lines, as needed by the DRP scheme. 
Despite some smearing of the wave front, the initial pulse, its propagation and scattering off the cylinder, 
were simulated fairly well (fig. 1). However, some oscillations inside the domain and spurious reflections 
off the boundaries started to emerge. 

Then, the grid topology was changed to an O-grid with a radius of 10.5D (cases 2 and 3). The grid 
had | -lines as concentric circles, with the first and last circles being respectively the cylinder and the outer 
boundary, and the 77 -lines emanated radially from the cylinder to the outer boundary (fig. 2). The time 
step At = 5.0E-3 was less than one-half of the accuracy limit for eq. (7) as applied to a linear wave 
equation, but about five times that needed for eq. ( 6 ). Although, the results with the 401 | -lines appeared 
adequate (case 2), doubling these lines rendered a truly symmetric initial pulse (case 3). On an SGI R10K 
computer in a time-shared mode, case 3 required 24 megabytes of run-time memory and 16 hours of CPU 
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processing. The unit processing time was computed to be about 0.2 ms/At/node. At three locations, given 
in their cylindrical coordinates, A(5D, ji/ 2), B(5D, 37i/4), C(5 d, tc), the pressure history was recorded from 
t = 6 to t = 10 at intervals of 0.01 (fig. 3). The computed results matched the analytical solution very 

well. The peak reached these points at about t ~ 6.3, 8.2 and 9.0, respectively, and it appeared slightly 
attenuating (0.06, 0.052, 0.048, respectively) due to the scattering. The peak-trough pair was followed 
by another set with lower amplitudes. All the waves were crisply simulated with virtually no numerical 
reflections from the boundaries (fig. 4). This very feature, i.e. the success of the boundary conditions, 
appeared to be pivotal for this problem. 

Periodic source: a limiting-cycle problem 


In the preparatory cases of 4 and 5, the equations were solved on an H-grid conforming to the wall shape: 
a flat plate in the former (fig. 5), and a circular bump on a flat plate in the latter case (fig. 6). The objective 
was to check the implementation of the boundary conditions and the suitability of the grid. The 
interference patterns from the cascades of incident and reflected waves reached a periodic state (limiting 
cycle), after some transient time, as could also be observed by the wall pressure (figs. 5 and 6). Note that 
the source amplitude was 1% in case 4. 

In cases 6 and 7, the scattering off - a cylinder was simulated on O-grids and the directivity, 


(16) 


D(6) = 




n 2 tTj Hj 


was computed at r=R and from time step nj to «2- In case 6, the solution was obtained on a 801x181 O- 
grid with 10.5D radius. This resulted in 20 PPW radially, and PPW circumferentially were: 28.6 on 
cylinder, 2.86 at r =5.0, and 1.36 at the outer boundary. However, from eq. (5), the theory required 4.5 
PPW, which was satisfied for points with r < 3.18. Consequently, despite the periodic response attained, 
e.g. on the cylinder (fig 7), neither the computed directivity pattern nor its amplitudes at i?=l0.5D were 
satisfactory. 


Hence, another O-grid was generated with 361x321 points and a radius of 8.5D (case 7). This resulted 
in 10 PPW radially, and PPW circumferentially were: 57.3 on cylinder, 5.73 at r =5.0, and 3.37 at the outer 
boundary. The theoretically required 4.5 PPW was satisfied for points with r < 6.36. (Practice may prove 
the safer requirement to be 8 PPW, which was satisfied for points with r < 3.58.) Initially, At was 2.5E-3, 
but after t=15, it was reduced to 1.25E-3, which was one-fourth of the accuracy limit for eq. (7) as applied 
to a linear wave equation, but 2.5 times that needed for eq. (6). Atr = R=5 and 6 from 0 to n/2, the 
periodic response was detected after 100 periods of source excitation, then the results were recorded at r 
=R =5 and 9 from 0 to 7t at 0.5-deg intervals. The computed directivity is presented in fig. 8. Although, 
the directivity had not yet attained the limiting cycle values at t =26.25, the number of peaks matched 
analytical values well. Since the results were relatively better for 9 < n/2, and for 9 > n/2, they improved 
with the elapsed time, it was deemed that all the transients had not yet left the domain. Also due to the 
marginal PPW circumferentially at r =5 and the uncertainty about the sufficient artificial viscosity to be used 
(Re in eq. (9) was set to 1.0E4), some oscillations were detected. This computed scattering pattern is also 
depicted via its pressure contours at two instants (fig. 9). Finally, on an SGI R10K computer in a time- 
shared mode, case 7 required 19.2 megabytes of run-time memory and 100 hours of CPU processing. 
Conceivably, the elapsed time for the scattering shown should have been doubled, which, naturally, 
would' have required twice as much computing time. 


Conclusions 


By and large, the present simulations of the propagation of acoustic waves, their reflections and scattering, 
in particular, the initial-value problem with the acoustic pulse, were successful. Two necessary building 
blocks to success, once a suitable CAA scheme was selected, were the correct boundary conditions, and an 
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adequate and efficient grid. Notwithstanding the imperfectly orthogonal grids and the required 
transformation metrics, employing the body-fitted coordinates allowed a straight forward implementation 
of the boundary conditions. The role of the grid became more accentuated in the periodic source case. 

The spread of the source (b in eq. (15)) and the intervals that the directivity was requested (1-deg) proved 
narrow enough to necessitate too fine a grid resolution, which superseded the benefits of a low PPW 
scheme. A better deployment of the grid points, such as, some sort of domain decomposition, could 
reduce the required computational resources. Also, since it took longer for the periodic behavior to be 
established at Q=n, it would have been less resource taxing to request the data up to, say, B=3nlA. Further, 
the definition of directivity included the /--><*>, leading one to place the outer boundary as far away as 
possible; however, the benchmark analytical solution was integrated virtually at any r -R value. Finally, a 
parametric study of the required amount of artificial dissipation proved to be another prerequisite. 
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Important characteristics of acoustic wave propagation are encoded in their disper- 
sion relations. Hence , a computational algorithm, which attempts to preserve these 
relations , was investigated. Considering the linearized , 2-D Euler equations, simula- 
tions were performed to validate this scheme and its boundary conditions . The results 
were found to agree favorably with the exact solutions . The boundary conditions 
were transparent to the outgoing waves, except when the disturbance source was 
close to a corner boundary. The time-domain data generated by such computations 
were often intractable until their spectra was analyzed. For this purpose, the relative 
merits of three spectral analysis methods were considered. For simple, periodic waves 
with steep-sloped spectra , the periodogram method produced better estimates than 
the Blackman-Tukey method, and the Hanning window was more effective when 
used with the former. For chaotic waves, however, the weighted-overlapped-segment - 
averaging and Blackman-Tukey methods were better than the periodogram method. 
Therefore , it was observed that the spectral representation of time-domain data was 
significantly dependent on the particular method employed. 


1 Introduction 

Computational aeroacoustics (CAA) may be defined as the 
employment of numerical techniques for the direct calculation 
of all aspects of aerodynamic sound generation and propagation 
starting from the time-dependent governing equations (Hardin, 
1993). Most computational fluid dynamics (CFD) schemes, 
however, are not adequately accurate for solving the aeroacous- 
tics problems (Lighthill, 1992). Their amplitudes are often or- 
ders of magnitude smaller, and yet the frequencies are orders 
of magnitude larger than the flow field variations generating the 
sound. Further, high-fidelity is paramount for the resolution of 
acoustic problems; but a consistent, stable, and convergent, 
high-order scheme is not necessarily dispersion-relation pre- 
serving and thus does not necessarily guarantee a good quality 
numerical wave solution for an acoustic problem. Hence, among 
the requirements that should be placed on a CAA algorithm are 
the minimal dispersion and dissipation features. Demonstrated 
in the present paper is an independent evaluation of a dispersion- 
relation-preserving (DRP) scheme, first introduced by Tam and 
Webb (1993) based on the group velocity considerations, and 
a compatible set of radiation and outflow boundary conditions. 

In the recent past, most commonly used computational meth- 
ods have been of the acoustic analogy type. With such a method, 
the aeroacoustic problem is reduced to a nonhomogeneous wave 
equation for the noise, with its right side equal to a distribution 
of acoustic sources of strength related directly to the characteris- 
tics of the flow (Farassat and Brentner, 1988). Another method 
is the Kirchhoff approach, where exploiting the similarities be- 
tween the aeroacoustic and the electrodynamics equations, the 
solution is obtained by integrating the wave equation external 
to some real or imaginary surface on which the relevant acoustic 
data is known (Lyrintzis, 1994). Perturbation methods for the 
linear potential equations have also been applied for the aeroa- 
coustic problems. Since a quiescent or a uniform flow is gener- 
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ally assumed, the potential methods are similar to the acoustic 
analogy equations except that the dependent variable in the 
wave /convective- wave operator is the perturbation velocity po- 
tential (Atassi et al., 1990). 

Perturbation Euler techniques have recently received atten- 
tion from the aeroacoustics (Hariharan and Bayliss, 1985) as 
well as the unsteady-flow communities (Meadows et al., 1993). 
The direct simulations of acoustic wave propagation have also 
been tried by solving the full Navier-Stokes equations. For ex- 
ample, Ridder and Beddini (1992) simulated the radiation of 
sound from resonance tubes, and Baysal et al. (1994) investi- 
gated two devices to suppress the high tones generated by a 
high-speed cavity flow. 

What is almost always common to these CAA solutions is 
the often intractable amount of time-dependent data generated. 
One way of reducing this data to understand the underlying 
physical phenomena is analyzing their spectrum, i.e., spectral 
analysis. A variety of methods are available to perform this 
task (Hardin, 1986). The spectra obtained from such analyses, 
however, may depend on the particular method that has been 
used. 

Therefore, considering this numerical data reduction as an 
integral part of CAA, three spectral analysis methods were in- 
cluded in the present study: Blackman-Tukey (1959), periodo- 
gram (Bartlett, 1950), and weighted-overlapped-segment-aver- 
aging (Welch, 1967) methods. The former two of these methods 
were implemented using box-car and Hanning windows. The 
latter method was compared for different amounts of data-seg- 
ment overlaps. Finally, these methods were applied not only 
for periodic but also for chaotic process data. For brevity, how- 
ever, some details reported by Vanel ( 1994) were not included 
herein. 


2 DRP Scheme 

The two-dimensional, compressible Euler equations were 
considered to simulate the propagation of waves in a uniform 
mean flow. By superimposing small perturbations on a mean 
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flow field, then by neglecting the higher order terms in perturbed 
quantities, these equations were linearized: 


dU dE dF 

s 7 + a 0 ^ = o 


(2.1. a) 


where the vector of unknowns ( U ) and flux vectors ( E , F) 
were 



~ p " 


PoU + pu Q 

u = 
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E = 

UqII + pi p 0 


u 


UqV 


„ P _ 


_ u 0 p + yp Q u _ 


Pou 



7PoV 


The density (p), pressure (p), and velocity («, v) of the per- 
turbed quantities are denoted without a subscript, but those of 
the mean flow are demarcated using the subscript 0. These 
variables were normalized using the following scales for length, 
velocity, time, density, and pressure, respectively: A* (mesh 
step, size), a 0 (speed of sound), A x/a 0t p 0 , and p^a\. 

Equations (2.1) support the acoustic, the entropy as well as 
the vorticity waves. The propagation characteristics of these 
waves (dispersion, dissipation, group and phase velocities, and 
isotropy or anisotropy ) are encoded in their dispersion relations, 
which relate the angular frequency of the waves (u>) to the 
wave numbers of the space variables (a and /?). 

In order to capture the correct wave propagation characteris- 
tics, the dispersion relation of the finite-difference scheme 
should match as closely as possible the dispersion relation of 
the partial differential equations (PDE). This is equivalent to 
requiring that the effective wave number (a) and angular fre- 
quency (ut) of the numerical scheme must be close approxima- 
tions to those of the PDE system for a large range of resolution. 
Such a scheme, therefore, is called a dispersion-relation pre- 
serving (DRP) scheme. In the remainder of this section, the 
way in which a DRP scheme was accomplished for the present 
computations is summarized. 

The spatial derivatives of the linearized Euler equations (2.1) 
were discretized using a central seven-point stencil. For exam- 
ple, the derivative df/dx was approximated at the It h node of 
i uniform grid as 


(df\ 1 2, 

l A Ax ajf,+j (2 - 2) 

expanding the right side of (2.2) in a Taylor series, five of the 
coefficients were determined by equating coefficients of the 
ame powers of A x (standard way), leaving two coefficients, 
ay a x and i, to be determined. Since for a central difference 
i = the resulting single parameter was determined by 

squiring that the numerical scheme and its corresponding origi- 
al set of governing equations would have the same dispersion 
elation. This process starts with the Fourier transform of the 
nite difference approximation and the partial derivative term 
right and left sides of (2.2), respectively). The result was 
te following relation between the numerical and exact wave 
umbers: 


“ = f- I a je' aiAx (2.3) 

j~~3 

ote that, aAx is a periodic function of aAx with a 2?r period, 
hen, a x was determined so that the integrated error £, between 
Ax and aAx was minimum: 
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E x = I |aAx - aAx\ 2 d(aAx) (2.4) 

J-rj 

where \aAx\ < rj . The integration limit 77 is determined de- 
pending on the deviation of Eq. (2.3) from the 45-deg line 
(for aAx > 1.45) and the shortest wavelength desired to be 
simulated. (Note that if Real(aAx) = (aAx), the scheme is 
nondispersive. and if Imaginary (5 Ax) = 0, the scheme is non- 
dissipative.) Here this length was chosen to be the minimum 
possible. 4.5 Ax, hence q = <t/ 2. (The choice of 77 = 1.1, 
which corresponds to the shortest wave length of 7 Ax, was 
later shown by Tam and Shen (1993) to produce better results 
for some waves.) This is a remarkable improvement over a 
standard, non-optimized. 6th-order scheme which can only re- 
solve wave lengths longer than 10 Ax. 

As for the time derivative OUldt, it was approximated by a 
four-level finite difference, which, in a standard sense, could 
be up to 3rd-order accurate: 

£/"' - U m =» A/ I J (2.5) 


where n is the time index. For the present scheme, however, 
the coefficients b,. b 2 and b, were determined in terms of b Q , 
so that when both sides were expanded in Taylor series, Eq. 
( 2.5 ) was satisfied to order ( At ) 2 . Taking the Laplace transform 
of the left and right sides of Eq. (2.5), but with t defined as a 
continuous variable, the effective angular frequency Q of the 
time discretization was found to be 


- 1 ) 

At 2 bje' JuA ‘ 
j- o 


( 2 . 6 ) 


The free parameter b„ was then chosen to minimize the inte- 
grated error E, between id At and uAt: 


= 0 


<r[Re(u)Af - ocAf)] 2 1 
+ (1 - ff)[Im(oAt - wAr)] 2 J 


d(uiAt) (2.7) 


where a is a weight that allows adjusting the compromise 
wished to impose on the optimization process between the dis- 
persion (real part) and dissipation characteristics (imaginary 
part). The present long-term time computations were performed 
with a = 0.36, which slightly biased the scheme in favor of a 
better dissipation character. 

After discretizing all the terms in Eq. (2.1), the resulting 
DRP scheme for the interior points was as follows: 

3 

U1*n = Uh« + A ‘ ^ b J K (2.8a) 

j = 0 

where 

Kl„ = - 4- E a,E%,„ - 4- S ajFl^j (2.8 b) 

Z\X^__i y y =-3 

/ and m are the spatial indices and n indicates the time level. 
The numerically stable maximum time step At was calculated 
from the Courant-Friedrichs-Lewy relation to be 


At 


CFL Ax 

1.75 [M + {1 + (Ax/ Ay ) 2 ) l/2 ] a 0 


(2.9) 


where M denotes the mean flow Mach number and CFL = 0.4. 
However, after analyzing the numerical damping of the time 
integration scheme and with cr — 0.36, the CFL value was set 
to the more stringent value ot 0. 19. 

Moving from the interior of the domain to the boundaries, 
they should be transparent to the disturbances reaching them to 
avoid any degradation of the numerical solution. For a right- 
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m oving uniform mean flow as represented by the linearized 
Euler equations, and when all the disturbances (acoustic, en- 
tropy and vorticity ) are generated in the interior of the domain, 
only acoustic waves reach the upper, lower and upstream bound- 
aries; but, in addition to the acoustic waves, entropy and vortic- 
ity waves can reach the downstream boundary. Therefore, fol- 
lowing Tam and Webb (1993) and from the asymptotic solu- 
tions of the finite difference form of Eq. (2.1 ), a set of radiation 
boundary conditions for all but the downstream boundary (in 
nondimensional form), 


dt dx dy 


(2.10a) 


where 


As J = 

r r 

■ c -£ 

r = (x 1 + y 2 ) U2 /Ax, V = 



(2.106) 

1/2 

(2.10c) 


and outflow boundary conditions for the downstream boundary, 


dU' 

dt 


+ M 


du ' 
dx 


s 


(2.11 a) 



x/ Ax 


<■) 


where 

(2.116) 


were obtained. In Eq. (2.1 1 ), M denotes the Mach number and 
U ' contains only the first three components of U in Eq. (2.16). 
The pressure was obtained as in Eq. (2.10). The points in the 
boundary regions were calculated by discretizing these equa- 
tions with the DRP approach: four-level time differencing and 
either a fully-backward or a fully-forward seven-point differ- 
ences spatially. 

The formal accuracy of the scheme was 4th-order in space 
and 2nd-order in time. The computational domain contained 
200 X 200 equally-spaced grid points, surrounded by three rows 
of grid points to form each of the boundary regions. 


5 = 


fy + M §p _dp _dp 
dt dx ’ dx ’ dy 


3 Wave Propagation Simulation 

The DRP scheme and its outflow and radiation boundary 
conditions were evaluated by considering three wave propaga- 
tion cases: (i) a single acoustic pulse, (ii) three successive 
acoustic pulses, (iii) two simultaneous acoustic pulses. Each of 
these pulses were introduced into a uniform mean flow (left to 
right) with a Mach number of 0.5. 

Each acoustic pulse was generated by setting u = v = 0 and 
imposing an initial Gaussian distribution for the pressure and 
density: 


p = p = 0.01 exp 


In 2 
(3Ax) 2 


(x 2 + y 2 ) 


(3.1) 


The computed pressure contours after 500 and 2000 time 
steps for the single acoustic pulse case are shown in Fig. 1. 
The results verified the expected propagation pattern: the radius 
of the acoustic wave expanded in time while its center was 
being entrained downstream with the mean flow. As shown in 
Fig. 16, the waves exited from the downstream, top, and bottom 
boundaries without any numerical reflections. Presented in Fig. 
2 is a comparison between the exact (Tam and Webb, 1993) 
and the numerical pressure wave forms at 500 and 2000 time 
steps. The numerical wave matches both the amplitude and the 
propagation speed of the exact wave. 



(b) 

Fig. 1 Selected pressure contours of a single acoustic pulse after (a) 
500 and (b) 2000 time steps. Intensity relative to initial peak disturbance 
pressure. 


To further test the effectiveness of the scheme and its bound- 
ary conditions, three successive acoustic pulses were generated 
in the center of the computational domain at t = 0, 600 A f, 
and 1200 Ar. Displayed in Fig. 3 are the pressure contours at 
1500 A t. Again, the wave centers were being drifted in the 
direction of the mean flow while their radii expanded continu- 
ously in time. The computation was carried out for a longer 
period of time only to observe that all three pulses left the 
computational domain, in succession, without any noticeable 
reflections. 

In the third wave propagation case, two acoustic pulses were 
generated simultaneously: one in the center of the computational 
domain and the other 50 Ax upstream of the outflow boundary 
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Fig. 2 Pressure wave form of a single acoustic pulse at the center. 


(as compared to the pulses generated 100 Ax away in the prior 
cases). The pressure contours at instants 500 At and 2000 At 
are plotted in Fig. 4. It was observed that the two pulses inter- 
sected as they propagated radially and drifted downstream. The 
interaction of the waves was crisply simulated. Small reflections 
started appearing (Fig. 4b) when both wave systems were exit- 
ing the outflow boundary. The primary suspect was the imple- 
mentation of the boundary conditions at the boundary comers. 
In fact, when the case was repeated on a four times larger 
domain (400 X 400 grid), now clear from the comers, the 
waves exited without any reflections (Fig. 5). Enlarging the 
domain, however, also placed the second source 150 A x up- 
stream of the outflow boundary. And, since the boundary condi- 
tions (Eqs. (2.10)-(2.11)) were developed from the asymp- 
totic solutions, their degradation for sources close to a boundary 
might not be entirely unexpected. 


4 Spectral Analysis Methods 

Two traditional methods for the spectral analysis of a station- 
ary process are the Blackman-Tukey (BT) and the periodogram 
(PM) methods. In the BT method, first the autocorrelation of 
the record of N data values was computed: 

1 N~r ~ i 

R(r) = -7 £ x(nAt)x((n + r)Af), 

N ~ T n = 0 

T = 0 T max (4.1) 

where r max was the maximum lag value obtained as the recipro- 




x/Ax 

Fig. 3 Selected pressure contours of successive acoustic pulses after 
1500 time steps. Intensity relative to initial peak disturbance pressure. 



(b) 

Fig. 4 Selected pressure contours of multiple acoustic pulses after: (a) 
500 and (6) 2000 time steps. Intensity relative to initial peak disturbance 
pressure. 
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Fig. 5 Selected pressure contours of multiple acoustic pulses after 2000 
time steps computed on an enlarged (400 x 400) domain. Intensity rela- 
tive to initial peak disturbance pressure. 


cal of twice the frequency resolution. Then, the power spectrum 
estimate was obtained by Fourier transforming Eq. (4.1): 

1 

P(f t ) = At[R(0) + 2 X R(r) cos (7r£r/r max ) } (4.2) 


where 


where 


box(r) = 



for 0 < t < T 

for t < 0 or t > T 


(4.6) 


Then, the Fourier transform of the boxcar window was a func- 
tion of type sin x/x, instead of a Dirac function, as it would be 
in the ideal case of having an infinite value of T. In order to 
reduce the side-lobe leakage, a window with small side lobes 
had to be used instead of the boxcar window. For this purpose, 
the Hanning window is commonly used: 


0.5 


h(t) = { 



0<r< T 


(4.7) 


0 , 


t < 0 or t > T 


Unfortunately, the side effect of side-lobe leakage reduction 
was always a broadening in the main lobe of the window trans- 
form. As a result, the main lobe of the Hanning window trans- 
form was twice as large as that of the boxcar window transform. 
This implied a reduction in the resolution quality by a factor 
of two. 

When the data of a chaotic process were to be analyzed, the 
problem of variability in the spectrum estimates, due to the 
limited length of the data record, also needed to be taken into 
consideration. A source of statistical error had been the time 
history tapering operation occurring in the PM method. It essen- 
tially discarded the relevant information near the beginning and 
the end of each record. To alleviate this increase in variability, 
the weighted-overlapped-segment-averaging ( WOSA) was em- 
ployed. In this method, the data record x(t) was broken into n d 
blocks, each consisting of N data values (Baysal et al., 1994). 
The power spectrum was then obtained by first taking the Fou- 
rier transform of the data over each block. 


f k = ~ and * = 0, 1 N/2 (4.3) 

NAt 

In the PM method, however, the power spectrum was ob- 
tained by taking directly the Fourier transform of the time do- 
main record to be analyzed. For a sampled data record r(«Ar), 
where n = 0, . . . , N — 1 , the power spectrum was 

N - 1 / 

At X x(nAt) exp I -y27r 

n— 0 \ 

where j = V— 1 and/* was defined as in Eq. (4.3). The highest 
frequency that could be reproduced from the data sampled at 
equal intervals At was the Nyquist frequency, f Nn . 

Three main problems arising when employing these methods 
were taken into consideration: the resolution , the leakage, and 
the statistical variability of the spectral estimates. A key param- 
eter in obtaining a reliable and accurate spectrum was the fre- 
quency resolution Af : the smaller the Af was ( longer record 
length) the more accurate and detailed the spectrum estimate 
became. A good resolution implied that each peak of the spec- 
trum estimate appeared sharply and distinctively. The leakage 
manifested itself by smearing the spectral estimate, i.e., energy 
in the main lobe of a spectral response leaked into the sidelobes, 
distorting other spectral responses that were present. This was 
due to the implicit windowing of the data that occurred when 
processing with the Fourier transform, as explained briefly be- 
low. 

The finite Fourier transform of a finite data record jt(r) with 
length Tmay be viewed as the Fourier transform of an unlimited 
time record u(r) multiplied by a boxcar window denoted by a 
function box (r): 


P(fk) - 


NAt 


nk 

N 


(4.4) 


x{t) = box ( r ) t) ( / ) 


(4.5) 


Xi(f k ) = At X Xi(n) exp( -j2n 


(n - IX* - 1) 
N 


for 

i = 1, . . . , n d ; n = U . . . , W; k = 1 AT/2 (4.8) 

then, averaging over the resulting spectra over all the blocks 


P{fk) = —~ f 

n d NAt 

|X ; (/*)| 2 


11 

to 

N! 2 

(4.9) 


Then, the data in each block was tapered by the Hanning win- 
dow to suppress the side- lobe leakage. 

5 Spectral Analysis Results 

The selected spectral methods for the present study were 
tested for their relative performance in analyzing periodic as 
well as chaotic processes. Since the waves considered in §.3 
would have trivial spectra, the periodic data record was gener- 
ated by a, multiple-frequency cosine function: 

F(t) ~ cos ( 27 t/o?) + cos ( 2-nfxt ) + cos (27r/ 2 f) 

+ 0.5 cos (2 tt/ 3 /) + 0.1 cos (2n f 4 t) (5.2) 

where the / values in kHz were:/ 0 = 1,/ = 1.1, f 2 = 2, / 3 = 
0.5,/ 4 = 0.2. The data were sampled at the frequency resolution 
of 10 kHz. Equation (5.2) was a good benchmark problem to 
illustrate the concepts of leakage and resolution, and to test the 
BT and PM methods in estimating such steep-sloped spectra. 
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Fig. 6 Power spectral density for multiple-frequency cosine function 
as represented by the periodogram method and different windows with 
resolution A f = 80 Hz. 


0 12 3 4 

Frequency (kHz) 

Fig. 8 Power spectral density for multiple-frequency cosine function 
as represented by the periodogram method and different windows with 
resolution Af = 50 Hz. 


For each method, the spectra were computed in four different 
ways by varying the resolution (A / = 80 Hz and 50 Hz) and 
by varying the windowing (boxcar and Hanning). 

In the PM results that are shown in Fig. 6, note the presence 
of leakage, i.e., non-zero values of the spectrum at frequencies 
off the peaks. These were reduced with the use of the Hanning 
window. Also, when the Hanning window was used, the off- 
peak SPL values were much lower (less leakage) than when the 
boxcar window was used. Figure 7 displays the power spectra 
obtained using the BT method. The distortion effect caused by 
the leakage appeared to be more prevalent as compared with 
the PM method. The large side lobes of the boxcar window 
were suppressed with the use of the Hanning window. However, 
the Hanning tapering yielded negative values and less side-lobe 
suppression in this case, as compared to its use with the PM 
method. From Figs. 6 and 7, it was also observed that the PM 
method produced four of the five peaks (/ 0 ,/j ,/ 2 ,/j); whereas, 
the BT method could only produce three of them (/o,/ 2 ,/ 3 ). 

In an attempt to demonstrate the effect of the frequency reso- 
lution, the spectra, initially obtained with A / = 80 Hz, were 


repeated with A / = 50 Hz. Shown in Figs. 8 and 9 are the PM 
and BT spectra, respectively, and they should be contrasted 
with Figs. 6 and 7. Lower frequency resolution resulted in more 
accurate spectrum representation. The effect of window tapering 
on the spectrum resolution was observed in Figs. 8 and 9: the 
peaks at f 0 and /i were not as distinct when the Hanning window 
was used. 

The next set of results were obtained for a chaotic process. 
The data record was obtained from Baysal et al. (1994): a 
computational simulation of a turbulent, transonic flow past a 
two-dimensional cavity with a length-to-depth ratio of 4.5. The 
instantaneous-pressure history was for a specific location inside 
the cavity. The SPL spectra were obtained using not only the 
WOSA (Figs. 10-13) but also the PM (Fig. 14) and the BT 
(Fig. 15) methods. An important parameter for the WOSA 
method is the amount of overlapping. This aspect was studied 
by comparing two different overlap cases with no overlapping 
(Figs. 11-13). 

The spectra shown in Figs. 10 and 11 were obtained by 
partitioning the data into three blocks with 50% of overlapping 



Frequency (kHz) 

Fig. 7 Power spectral density for multiple-frequency cosine function 
as represented by Blackman-Tukey method and different windows with 
resolution A/ = 80 Hz. 


Frequency (kHz) 

Fig. 9 Power spectral density for multiple-frequency cosine function 
as represented by Blackman-Tukey method and different windows with 
resolution Af - 50 Hz. 
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Fig. 10 Sound pressure levels in a cavity as represented by WOSA Fig. 13 Sound pressure levels in a cavity as represented by WOSA 
method and boxcar window with 50 percent overlapping. method and Hanning window with 75 percent overlapping. 
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Fig. 11 Sound pressure levels in a cavity as represented by WOSA 
method and Hanning window with 50 percent overlapping. 
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Fig. 12 Sound pressure levels in a cavity as represented by WOSA 
method and Hanning window with no overlapping. 


between them. The spectral peaks appeared more clearly (espe- 
cially, for higher frequencies) when the Hanning window taper- 
ing was done. It was also noted, that the broadband SPL de- 
creased in intensity with increasing frequency, whereas it re- 



Frequency (kHz) 

Fig. 1 4 Sound pressure levels in a cavity as represented by periodogram 
method and Hanning window. 



Frequency (kHz) 

Fig. 15 Sound pressure levels in a cavity as represented by Blackman- 
Tukey method and Hanning window. 


mained roughly constant when a boxcar window was used. 
Comparing the spectra of 50% overlapping with no overlapping 
(Fig. 12), it was seen that the no-overlapping spectrum con- 
tained more local peaks. Also, since Fig. 1 1 displayed crisper 
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peaks and less variation, it was concluded that overlapping in- 
creased the stability. However, almost no distinguishable differ- 
ences appeared between the spectrum with 50% overlapping 
and the one obtained with 75% overlapping (Fig. 13). Since 
the amount of overlapping increased the number of Fourier 
transform operations, using more than 50% overlapping did not 
appear to be feasible. This was in agreement with Otnes and 
Enochson ( 1972) that 50% overlapping retrieves about 90% of 
the stability lost due to the tapering operation, hence it has been 
more commonly used. 

When the PM method with a Hanning window was used 
(Fig. 14), more local peaks were generated (as compared to 
Fig. 12), and only the ones located in the lower end of the 
spectrum were clearly distinguishable. However, the spectrum 
obtained by using the BT method with a Hanning window (Fig. 
15 ) compared very well with the ones obtained using the WOS A 
method and 50% overlapping. 

6 Conclusions 

The present results showed that the DRP scheme was reliable 
to yield accurate simulations of isotropic waves with wave- 
lengths as small as 4.5 Ajc. The radiation and outflow boundary 
conditions were very effective unless a disturbance source was 
placed close to a boundary comer. 

In order for the scheme to be of more practical interest, 
it was deemed necessary to investigate its extensions for: (i) 
nonuniform grids, (ii) multiblock grids, (iii) time integration 
requiring less computer storage, and (iv) for the nonlinear equa- 
tions with a controllable amount of numerical dissipation. Also, 
both from the efficiency and the accuracy points of view, replac- 
ing the present finite-difference (local) boundary conditions 
with integral (global) boundary conditions should prove to be 
profitable. 

For the spectral analysis of the time-domain data obtained 
from a periodic process with a steep-sloped spectrum, the PM 
method produced better estimates than the BT method on the 
basis of minimum leakage. The Hanning window, when used 
with the BT method, was found to yield less leakage suppression 
than when it was used with the PM method. 

For a chaotic process, on the other hand, the WOSA and BT 
methods were shown to yield stable results of comparable qual- 
ity. In the WOSA method, the extra computational cost of over- 
lapping more than 50% could not be justified. Although the 
methods were not methodically studied for their relative compu- 
tational efficiencies, merely based on the operation count, it 
was estimated that the WOSA method would be more economi- 
cal than the BT method. 

Finally, it was demonstrated that the very same time-depen- 
dent data, whether it be periodic or chaotic, could easily be 
interpreted discrepantly as for its quality and the underlying 


physical message, depending on the spectral analysis method 
employed. 
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Abstract 

Accurate computation of acoustic wave propagation may 
be more efficiently performed when their dispersion relations are 
considered. Consequently, computational algorithms which 
attempt to preserve these relations have been gaining popularity 
in recent years. In the present paper, the extensions to one such 
scheme are discussed. By solving the linearized, 2-D Euler and 
Navier-Stokes equations with such a method for the acoustic 
wave propagation, several issues were investigated. Among 
them were higher-order accuracy, choice of boundary conditions 
and differencing stencils, effects of viscosity, low-storage time 
integration, generalized curvilinear coordinates, periodic sources, 
their reflections and interference patterns from a flat wall and 
scattering from a circular cylinder. The results were found to be 
promising en route to the aeroacoustic simulations of realistic 
engineering problems.. 

Introduction 

Computational aeroacoustics (CAA) may be defined as the 
application of numerical techniques for the direct calculation of 
aerodynamic sound generation and propagation starting from the 
first principles. Most computational fluid dynamics (CFD) 
schemes, however, are not adequately accurate for solving the 
aeroacoustics problems (Lighthill, 1992). Their amplitudes are 
often orders of magnitude smaller, and yet the frequencies are 
orders of magnitude larger than the flow field variations 
generating the sound. Further, high-fidelity is paramount for the 
resolution of acoustic problems; but a consistent, stable, and 
convergent, high-order scheme is not necessarily dispersion- 
relation preserving and thus does not necessarily guarantee a 
good quality numerical wave solution for an acoustic problem. 
Hence,, among the requirements that should be placed on a CAA 
algorithm are the minimal dispersion and dissipation features 
(Tam and Webb, 1993). 

The direct simulations of acoustic wave propagation have 
been tried by solving the full Navier-Stokes equations. For 
example, Baysal et al. (1994) investigated two devices to 
suppress the high tones generated by a high-speed cavity flow. 


using an unsteady computational fluid dynamics. (CFD) method. 
The comparisons with experimental data were acceptable for 
engineering purposes. However, it was realized that the 
dissipative and dispersive characteristics of a typical second - 
order CFD method would preclude a long-term wave propagation 
simulation. Also, various studies (e.g. Mankbadi et al., 1993, 
Lyrintzis et al., 1995, and Hardin et al., 1995) have suggested that 
the direct simulations of the flow equations for the acoustic wave 
propagation using the higher-order CFD schemes could become 
prohibitively expensive, since the number of grid points per 
wavelength would be excessively high (ideally should not exceed 
ten). 

Therefore, a fourth-order accurate dispersion-relation - 
preserving (DRP) method was previously investigated for a 
variety of wave propagation problems (Vanel and Baysal, 1995). 
A number of observations and recommendations were made for 
future investigations to extend the scheme to solve some 
application problems. The present investigation started precisely 
with ‘this impetus. A selection of the issues explored, all for 
linear cases, are reported herein. These include higher-order 
accuracy (up to sixth-order), choice of boundary conditions and 
differencing stencils, effects of viscosity, low-storage time 
integration, generalized curvilinear coordinates, periodic sources, 
their reflections and interference patterns from a flat wall and 
scattering from a circular cylinder. The extensions for the 
nonlinear acoustics are deferred to another paper for brevity. 

Governing Equations 

The two-dimensional, compressible Navier-Stokes and 
Euler equations were considered in generalized curvilinear 
coordinates. In the absence of curved or irregularly shaped 
boundaries, their Cartesian expressions were preferred. By 
superimposing small perturbations on a mean flow field, then by 
neglecting the higher order terms in perturbed quantities, these 
equations were linearized to simulate the propagation of waves in 
a uniform mean flow.in x-direction. 

■rj— = -R(U) + S .where R(U) = ^+-^ (l.a) 

9t dt 8t| 
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The flux vectors ( E,F) were obtained through transformations. 


E = [4 x (E-E v )+^ y (F-F v )]/J. 

F = [ii x (E-E v ) + Ti y (F-F v )]/J (l.b) 

where the vector of unknowns ( U ) and the physical fluxes were 
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The diffusion terms in the energy equation were omitted due to 
their perceived secondary influence on most of the acoustic wave 
propagation phenomena. Also, the second coefficient of 
viscosity was neglected, effectively negating the rotational effect 
of a fluid element but leaving its deformation-rate in the shear 
stresses. 


x 


_ H» r duj , 

dxj dx, 


(l.e) 


The density (p), pressure (p), and velocity (u,v) of the perturbed 
quantities are denoted without a subscript, but those of the mean 
flow are demarcated using the subscript 0. These variables were 
then normalized using the following scales for length, velocity, 
time, density, and pressure, respectively: Ax (mesh step size), 
a 0 (speed of sound), Ax/a 0 p 0 , and p 0 aQ. Finally, S in eq. 
(1) denotes a possible acoustic source. 

Computational Method 

Equation (1) supports the acoustic, the entropy as well as 
the vorticity waves. The propagation characteristics of these 
waves (dispersion, dissipation, group and phase velocities, and 
isotropy or anisotropy) are formulated in their dispersion 
relations, which relate the angular frequency of the waves (co) to 
the wave numbers of the space variables (a). Therefore, in order 
to capture the correct wave propagation characteristics, the 
dispersion relation of the finite-difference scheme should match 
as closely as possible the dispersion relation of the partial 
differential equations (PDE). This is equivalent to requiring that 
the effective wave number (a) and effective angular frequency 
( to ) of the numerical scheme must be close approximations to 
those of the PDE system for a large range of resolution. Such a 
scheme, therefore, is called a dispersion-relation preserving 
(DRP) scheme. The way in which the baseline DRP scheme of 
Tam and Webb (1993) and Tam and Shen (1993) was 
accomplished for the present development is given by Vanel and 
Baysal (1995). 


In discretizing the spatial derivatives, the DRP was 
achieved by determining the coefficients (aj) from the Taylor 
series expansion as a one-parameter family; then the remaining 
coefficient was determined from the minimization of the 
discrepancy between the numerical (effective) wave number a 
and the exact wave number a , by integrating the square of their 
difference for a desired range (-£,£), where |aAx|^£. Note 
that, aAx is a periodic function of aAx with a 2k period, and if 
Real(a Ax) = (a Ax), the scheme is nondissipative, and if 
Imaginary(a Ax) = 0, the scheme is nondispersive. The 
integration limit £ is determined depending on the shortest 
wavelength desired to be simulated. For example, in the 4-th 
order version of the present schemes, two choices of £ were tried: 
7t/2 and 1.1, which corresponded to the minimum wavelengths of 
4.5 Ax and 7 Ax, respectively. This is a remarkable improvement 
over a standard sixth-order scheme which can only resolve wave 
lengths longer than 10 Ax. 

The time integration of eq. (1) was performed in two 
different ways. In the first approach, a four-point finite 
difference, which in a standard sense could be up to third-order 
accurate, was derived from the Taylor series as a one-parameter 
family. The remaining coefficient (bj) was determined, as in the 
spatial coefficients, by minimizing the discrepancy between the 
effective and the exact dispersion relations. After discretizing all 
the terms in eq. (1), the resulting d( At 2 , Ax N+M “ 2 ) -accurate 
DRP scheme was as follows: 

U"m=U?. m +At ibjR^-j (2.a) 

j“0 

I M 1 M 

whef e I a jF"m+j (2-b) 

Aq j=-N J All j=-n j 

t and m are the spatial indices and n indicates the time level. 
For N=M, difference eq. (2.) is central, for N=0, it is fully 
forward, and for M=0, it is fully backward. All the interior cells 
were computed using central differences. However, since these 
high-order stencils require multiple layers of boundary cells, all 
combinations between a central and a fully-one-sided difference 
need also be derived. Only then, it would be possible to always 
utilize the information from the nearest possible points for better 
accuracy. In the present computations a fourth-order scheme and 
a sixth-order scheme were derived, requiring 7-point stencil (N 
takes values from 6 to 0, M takes values from 0 to 6, and 
N+A/=<5) and 9-point stencil {N takes values from 8 to 0, M takes 
values from 0 to 8, and N+M=8) , respectively. 

The numerically stable maximum time step At was 
calculated from the Courant-Friedrichs-Lewy relation. For 
example, for the fourth-order scheme in Cartesian coordinates, 
the stable CFL number was found to be 0.4. However, after 
analyzing the numerical damping of the time integration scheme, 
the CFL value was set to the more stringent value of 0. 1 9. 

Since the above time integration scheme required the 
storage of four time levels, a relatively lower storage alternative, 
such as the Runge-Kutta scheme, was considered. The classical 
Runge-Kutta schemes, however, are intrinsically dissipative and 
dispersive. Recently, a class of low-dissipation and low- 
dispersion Runge-Kutta schemes have been developed by Hu et 
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al. (1994). Its development was similar to that of eq. (2), such 
that the dissipation and dispersion errors were minimized for all 
the frequencies resolvable by the numerical discretization. Most 
importantly, these schemes can be implemented with low-storage 
requirements. The resulting scheme had the spatial integration 
identical to eq. (2b), but the time integration was replaced by the 
following: 


where the Mach-number-related operator and the source term 
were. 


d 3 

M=(u 0 ^ x +v 0 ^ y )— +(u 0 T 1 x + v 0 Tl y )— (6b) 


Q = 


dP , \A 

at 


dp j_ dp> 


dp, n dp. 


-|T 


+ M(p), 


(6c) 


U (0) =U n (3a) 

U (i) =U (0) -3 i AtR (i_l) , i = 1,2,.., p (3b) 

U" +1 =U (P) (3c) 

The indices n, p and i indicate the time level, and the order and 
the stage of the Runge-Kutta method, respectively. As for the 
coefficients, and the other coefficients pj were determined 
from 


were obtained. In eq. (6), U' contains only the first three 

components of U in Eq. (1c). The pressure was obtained as in 
eq. (5). 

For the inviscid calculations on a solid wall , the 
impermeability condition requires that the normal contravariant 
velocity be zero; and for the viscous computations, the no-slip 
condition requires that the tangential contravariant velocity also 
be set to zero. 


Ci=nPp-k+ 2 * i = 2,...,p 

k=2 


(4) 


The coefficients q were computed by considering the 
amplification factor of the Runge-Kutta scheme, then minimizing 
the dispersion-relation error. The time steps to be used were 
determined from the stability as well as the accuracy limits. In 
the present study, a five-stage Runge-Kutta (p=5) was used, 
which required two levels of storage and it was at least second- 
order accurate. When it was used with 7-point spatial stencil, the 
CFL limit from the stability was found to be 3.05, but it was only 
1.16 from the accuracy limit. Since, however, this still was larger 
than the CFL limit of the DRP time integration (0.19), this 
method was also more efficient in processing time. 

Boundary Conditions 

The boundaries should be transparent to the acoustic 
disturbances reaching them to avoid any degradation of the 
numerical solution. For a right-moving uniform mean flow as 
represented by the linearized equations, and when all the 
disturbances (acoustic, entropy and vorticity) are generated in the 
interior of the domain, only acoustic waves reach the upper, 
lower and upstream boundaries; but, in addition to the acoustic 
waves, entropy and vorticity waves can reach the downstream 
boundary. Therefore, following Tam and Webb (1993) and from 
the asymptotic solutions of the finite difference form of eq. (1), a 
set of radiation boundary conditions. 


au . au D au ~ n 

— +A— +B— +CU==0 

at a^ an 


where 


(5a) 


x^ x -**y£ v xn x +yTiv v 

A = V y , BsV — ■ y , C=— (5b) 

r r 2r 

r=— -Jx 2 +y 2 , V=— M+,|1-(M— ) 2 , (5c) 

Ax r V r 

and outflow boundary conditions for the downstream boundary. 


^+M(U') = Q. 


(6a) 


V = T| X U +T}yV = 0. , U = ^ x u4*^ y v = 0 (7) 

When above equations were used in the n -momentum and £- 
momentum equations, the wall values of pressure and shear stress 
were obtained, respectively. The DRP scheme coefficients for all 
the boundary conditions (Tam and Dong, 1994) were derived by 
an analogous method to that of the boundary region cells. 

Results 

The present schemes and their boundary conditions were 
evaluated by considering five wave propagation cases: (i) Single 
acoustic pulse; (ii) Two simultaneous acoustic pulses; (iii) 
Acoustic pulse near a flat wall; (iv) Acoustic pulse near a 
circular cylinder; (v) Periodic source near a flat wall; (vi) 
Periodic source near a flat wall with a circular bump. 

Each of the acoustic pulses, introduced into a uniform mean 
flow (left to right) with a Mach number of 0.5, was generated by 
setting u=v=0, and imposing an initial Gaussian distribution for 
the pressure and density: 


p = p exp {-In2 [ 


(x-x s ) 2 +(y-y, s ) 2 


( 8 ) 


For the cases (i)-(iv), p = 0.01,b = 3Ax, x s =y s =0. The 
pressure contours at 1 500 At for the single acoustic pulse case, 
computed by the fourth order inviscid DRP scheme (eq. 2), are 
shown in Fig. la. The results verified the expected propagation 
pattern and matched the exact solution (Vanel and Baysai, 1995): 
the radius of the acoustic wave expanded in time while its center 
was being entrained downstream with the mean flow. The waves 
exited from boundaries without any numerical reflections. 
Presented in Fig. lb is the same case computed using the 
optimized Runge-Kutta time integration (eq. 3) after only 150 
time steps. The numerical waves matched in both the amplitude 
and the propagation speed. However, the latter required about 
one-half the storage memory and one-sixth the processing time of 
the former. 

In the second case, two acoustic pulses were generated 
simultaneously: one at half-span and the other at three-quarter 
span of the computational domain. Presented in Fig. 2a are the 
pressure contours computed by the fourth-order DRP scheme 


3 



Proc. ASME Fluids Engineering Division Summer Meeting, 1996 

FED-Vol. 238, pp. 503-510. 
Kaushik and Baysal 


with e=7t/2 and without a mean flow (M=0). It was observed that 
the two pulses intersected as they propagated radially. The 
interaction of the waves was crisply simulated. Then, the case 
was repeated with a mean flow of M=0.5 and £=1.1 (Fig. 2b). 
With the mean flow, the waves started to drift downstream. 
Then, the case was repeated using the sixth-order scheme (Figs. 
2c, 2d). In all of the above cases minor reflections from the 
comers were observed. This was cured when the simulation was 
repeated on a four times larger domain (Figs. 2e, 2f). 

To test the wall boundary condition, the lower boundary 
was replaced by a solid flat wall. The reflection of the pulse off 
the wall and its interference with the incident pulse were 
computed once by solving the Euler equations (Fig. 3), then by 
solving the Navier-Stokes equations (Fig. 4). As expected, the 
wave strength was attenuated faster in the latter due to the 
physical diffusion process. 

Then, the case of circular cylinder, with the pulse at (4, 0) 
with p = 1.0,b = 0.2, was simulated on a relatively coarse but 
conforming grid (251x101). The motivation was the physical 
problem of the sound field generated by a propeller scattered off 
by the fuselage of an aircraft. The pressure loading on the 
fuselage was an input to the interior noise problem. Here, the 
fuselage was idealized as a circular cylinder and the noise source 
as a line source, hence a 2-D problem. Four instants from the 
animation of the acoustic scattering are presented in Fig. 5. 

In the periodic acoustic source cases, the source was 
generated with the fourth term in the source vector of eq. (1), 

S 4 = s exp{-ln2 + (y&i)-y») 2 ft sin(cot) (9) 

b 

where b=3, and (x s ,y s ,)=(0,20). The medium was inviscid air at 
M=0. The interference pattern, shown in Fig, 6 is for the source 
with s = 0.001 near the flat wall, and after the computations 
reached the limiting cycle. Then, the periodic source with 
s » 0.01 was placed above a circular bump on a flat plate at (0, 
20). The interference pattern is presented in Fig. 7. 

By and large, the reflections and scattering 
computations were successful. By employing the body-fitted 
coordinates, solid wall boundary conditions were imposed 
properly, hence the wall region was properly computed. On the 
other hand, grid cells not being perfectly orthogonal, and the 
existence of the transformation metrics in the equations, caused 
minor degradation, as compared to the Cartesian cases. 

Conclusions 

Computational schemes which preserve the dispersion 
relation of the fundamental equations were investigated. In 
developing such schemes, ultimately for realistic and complex 
aeroacoustics problems, several necessary steps have been taken. 
The linearized Euler and Navier-Stokes equations were solved. 
Solid wall boundary conditions and differences utilizing the 
nearest-point information in the boundary regions were 
demonstrated. The spatial scheme was extended to the sixtth 
order accuracy, and with the optimized Runge-Kutta time- 
integration, the efficiency of the method was improved. For the 
curved wail boundary conditions, the scheme was extended for 
the generalized curvilinear coordinates. The schemes were 
successfully demonstrated for several acoustic pulse cases and 
several acoustic scattering cases. 


In order for the scheme to be of more practical interest, it 
was deemed necessary to investigate its further extensions for the 
nonlinear equations, with a controllable amount of numerical 
dissipation, and multiblock nonuniform grids. Also, both from 
the efficiency and the accuracy points of view, replacing the 
present finite-difference (local) boundary conditions with integral 
(global) boundary conditions should prove to be profitable. 

References 

Baysal, O., Yen, G.W., and Fouladi, K., 1992, “Navier-Stokes 
Computations of Cavity Aeroacoustics With Suppression 
Devices,” Proceedings of DGLRJAIAA 1 4th Aeroacoustics 
Conference, Vol. 2, pp. 940-948, Aachen, Germany. Also, 
Journal of Vibration and Acoustics, Vol. 116, No. 1, pp. 105- 
112. 

Hardin, J.C., Ristorcelii, J.R., Tam, C.K.W., (Editors), 1995, 
ICASE/LaRC Workshop on Benchmark, .Problems in 
C omputational AcroacQ .u.s fc .$, NASA Conference Publication 
3300. 

Hu, F.Q., Hussaini, M.Y., and Manthey, J., 1994, "Low- 
dissipation and Low-dispersion Runge-Kutta Schemes for 
Computational Acoustics," ICASE Report 94-102, Hampton, 
VA. Also to appear in Journal of Computational Physics , 
March 1996. 

Lighthill, J., (1992), "Report on the Final Panel Discussion on 
Computational Aeroacoustics," ICASE Report No. 92-53, 
NASA Langley Research Center, Hampton, VA. 

Lyrintzis, A. S., Mankbadi, R.R., Baysal, O., Ikegawa, M., 
(Editors), 1995, Computational Aeroacoustics, , FED-Vol. 219, 
ASME, New- York, NY. 

Mankbadi, R.R., Lyrintzis, A. S., Baysal, O., Povinelli, L. A., 
Hussaini, M. Y., (Editors), 1993, Computational Aero- and 
Hydro-acoustics . FED-Vol. 147, ASME, New-York, NY. 

Tam, K. W., and Webb, J. C., 1993, "Dispersion-Relation- 
Preserving Finite Difference Schemes for Computational 
Acoustics," Journal of Computational Physics , Vol. 107, pp. 
262-283. 

Tam, C.K.W., and Shen, H., 1993, "Direct Computation of 
Nonlinear Acoustic Pulses using Higher-Order Finite 
Difference Schemes," AIAA Paper 93-4325, 15th 
Aeroacoustics Conference, Long Beach, CA. 

Tam, C.K.W., and Dong, Z., 1994 "Wall Boundary Conditions 
for High-Order Finite Difference Schemes in Computational 
Aeroacoustics," AIAA Paper 94-0457, 32nd Aerospace 
Sciences Meeting, Reno, NV. 

Vanel, F. O., and Baysal, O., 1995, “Investigation of Dispersion- 
Relation-Preserving Scheme and Spectral Analysis Methods 
for Acoustic Waves,” Paper no. 95-093, Proceedings of First 
CEAS/AIAA Aeroacoustics Conference . Munich, Germany, 
pp. 675-682. Also, to appear in [ASME] Journal of Vibration 
and Acoustics , 

Acknowledgment 

This work was supported by NASA Langley Research Center 
Grant NAG-1- 1653. The technical monitor was Dr. J.L. Thomas. 


4 



- 5 < 


Fig. 1 Pressure contours of an acoustic pulse (eq. 8) after 1500 At computed with: (a) four-level DRP time 

integration, (b) optimized five-stage Runge-Kutta. 



Fig. 2 Pressure contours of multiple acoustic pulses (eq. 8 at two locations), (a) M=0 case at 1000 At;(b) e=1.1 case 
at 1500 At; Sixth-order scheme at (c) 1500 At and (d) 2000 At; On 400x400 domain at (e) 1000 At and (f) 2000 At. 
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(c) 

Fig. 3 Reflection of an acoustic pulse (eq. 8) from a flat wall 
(a) 500 At, (c) 700 At, (d) 1000 At. 


(d) 

by solving Euler equations. Pressure contours at: 
(b) Velocity vectors at 500 At. 
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(a) (b) 

Fig. 6 Interference pattern of a periodic acoustic source (eq. 9) reflected from a flat wall: 
(a) pressure contours, (b) wall pressure values. 




(a) (b) 

Fig. 7 Interference pattern of a periodic acoustic source (eq. 10) from a bump on a flat wall: 
(a) pressure contours, (b) wail pressure values. 
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THE ACCURACY OF SHOCK CAPTURING IN TWO SPATIAL DIMENSIONS 


Mark H. Carpenter * Jay H. Casper * 


Abstract 

An assessment of the accuracy of shock captur- 
ing schemes is made for two-dimensional steady flow 
around a cylindrical projectile. Both a linear fourth- 
order method and a nonlinear third-order method 
are used in this study. It is shown, contrary to con- 
ventional wisdom, that captured two-dimensional 
shocks are asymptotically first-order, regardless of 
the design accuracy of the numerical method. The 
practical implications of this finding are discussed 
in the context of the efficacy of high-order numeri- 
cal methods for discontinuous flows. 

Introduction 

High-order schemes (schemes of second-order ac- 
curacy or higher) have been used effectively to cap- 
ture shocks for at least thirty years 1 . There is lit- 
tle doubt of their virtues when compared with first- 
order schemes. For this reason, second- and higher- 
order schemes have almost universally replaced first- 
order schemes throughout computational fluid dy- 
namics. Only recently, however, has a critical anal- 
ysis of the actual error of higher-order schemes been 
undertaken. In a previous work by Casper and Car- 
penter 2 it is shown that time dependent shocks can 
only be captured with first-order accuracy, regard- 
less of the design order of the numerical method. 
Figure 1 shows the pointwise error from a sound- 
shock interaction problem, using a fourth-order non- 
linear ENO algorithm to solve the time-dependent 
Euler equations. Figure 2 shows the same problem 
resolved using a sixth-order linear finite difference 
algorithm. Both algorithms converge at a first-order 
rate, downstream of the shock. This degradation in 
accuracy is caused by the inability of conventional 
numerical methods to pass information through a 
discontinuity. 
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Figure 1 . Pointwise error for the 1-D sound-shock 
interaction problem as obtained using the nonlinear 
ENO-4-3 numerical scheme . 

In a companion work these results are reaf- 
firmed and quantified for a broader class of prob- 
lems. Specifically, a sixth-order- accurate compact 
implicit finite difference scheme and a fourth-order 
accurate ENO scheme are used to investigate various 
discontinuous flows. The design order of accuracy 
is achieved in the smooth regions of a steady-state, 
quasi-one-dimensional Euler test case, as well as in 
the time-dependent Burgers’ equation. However, in 
an unsteady Euler sound-shock interaction problem, 
first-order results are obtained downstream of the 
shock. A discontinuous linear model problem identi- 
fies the cause of the first-order results, and quantifies 
the error as being predominantly a numerical phase 
shift resulting from information passing through the 
discontinuity. Post-processing the first-order data, 
increases the solution accuracy downstream of the 
discontinuity to second-order, although high-order 
is still not attainable. 

One can readily extend the negative results found 
in 1-D to multiple spatial dimensions. Specifi- 
cally, that time-dependent shocks are inherently 
first-order accurate in any number of spatial dimen- 
sions. The 1-D steady Euler equations exhibit de- 
sign accuracy away from discontinuities. There is, 
therefore, the possibility that the steady-state dis- 
continuous Euler equations will admit higher-order 
solutions in multiple dimensions. In this work we 
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Figure 2. Pointwise error for the 1-D sound-shock 
interaction problem as obtained using linear LIN-6-4 
numerical scheme . 


quantify the solution accuracy for the steady-state 
2-D Euler equations. We show that, unlike the 1-D 
case, the 2-D Euler equations are first-order accurate 
downstream of a shock at steady-state, regardless 
of the design accuracy of the numerical method. 

In section 2, we present the Euler Equations in 
two spatial dimensions. The equations are used to 
solve for the inviscid flow around a supersonic blunt 
body. In section 3, we describe the numerical meth- 
ods used on the test problems. Three different nu- 
merical methods are used to solve for the flow around 
a supersonic blunt-body. A Chebyshev bow-shock 
fitting algorithm is used to obtain an “exact solu- 
tion” to the problem. A linear fourth-order finite dif- 
ference algorithm, and a nonlinear third-order ENO 
numerical algorithm are used to capture the solution 
around the blunt-body. In section 4, we present a 
comparative study between the three methods. We 
show that all captured solutions are first-order on 
sufficiently fine meshes, regardless of the design ac- 
curacy of the spatial operator. In Section 5, we study 
the influence of Mach number, and the effects of de- 
sign accuracy on the first-order results. Section 6 
concludes the work. 

Governing Equations 

We focus in this work on the two dimensional Eu- 
ler equations. In the physical coordinates the equa- 
tions are: 


where 
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pu 


pv 

u= 

pu 
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pu 2 + P 
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pvu 


pv 


puv 


pv 2 + P 




_(pE + P)u_ 


~(pE 4- P)v _ 


The variables p ) u, v y P, E , are the density, x-velocity, 
y-velocity, pressure, and total specific energy, respec- 
tively. The equation of state is 

p = (7~ 1 )p[E- ^(u 2 + v 2 )] 

where 7 is the ratio of specific heats, which is as- 
sumed to have a constant value of 1.4. 

Numerical Methods 

We use three different numerical methods to study 
the 2-D blunt body problem. A Chebyshev shock fit- 
ting algorithm is used to generate an “exact” solu- 
tion around a supersonic cylindrical projectile. Both 
a linear fourth-order and a nonlinear third-order 
ENO shock capturing algorithm are then used to 
capture the solution around the cylindrical body. 
Accuracies are assessed by comparing the captured 
solutions with the exact solution. All three of these 
algorithms are fairly routine, and have been doc- 
umented extensively in the literature. We include 
only a brief description of each algorithm, noting 
the specific details necessary to maintain robustness 
on the blunt body problem. 

When dealing with negative results, one must be 
wary of the generality of the test problems, and of 
the numerical algorithms. The two capturing al- 
gorithms represent extremes from the spectrum of 
higher-order formulations presently in use. It is 
hoped that this diversity adds credibility to the neg- 
ative results obtained in this study. 

Chebyshev Shock Fitting 

A Chebyshev bow-shock fitting algorithm is used 
to determine a very accurate numerical solution to 
the blunt body problem. Figure 3 shows a numer- 
ical grid wrapping around the forebody of a two- 
dimensional circular cylinder. A Chebyshev spatial 
operator is used in both the radial and the circumfer- 
ential directions to resolve the inviscid flow around 
the cylinder. Note that the grid line furthest from 
the body coincides exactly with the curved bow- 
shock. 

The fitted solution is obtained by first mapping 
the Euler equations in space and time into the com- 
putational space. A linear mapping from (r,0,t) — > 
(£, 77, r) is used: 





Figure 3. Shock- fit Chebyshev spectral grid for the 
Mach = 2.5 blunt-body flow . 


Figure 4. Pressure profiles of the “exact” solution 
obtained with the Chebyshev shock fitting algorithm . 


V 

r 

where r by r 5y 6 m in > and 9 max , are the radius of the 
body, and shock, and the minimum and maximum 
values of 9 , respectively. The shock position always 
coincides with the line £ = 1, and the physical do- 
main changes continuously as the bow-shock moves 
to its steady state position. Because of symmetry 
along the body centerline, only half of the prob- 
lem is solved. The Chebyshev collocation technique 
is used in the computational domain to discretize 
the resulting set of equations. The equations are 
then marched in time until a steady state solution is 
reached. 

A crucial element of the procedure is the treat- 
ment of the boundary conditions, and in particular 
the movement of the shock. To move the shock, a 
shock velocity equation is derived by differentiating 
the Rankine-Hugoniot relations at the shock. The 
resulting differential equations for the position and 
velocity of the shock are then advanced in time with 
a five stage Runge-Kutta scheme along with the rest 
of the discrete equations. Symmetry conditions are 
imposed along the center line, and in viscid wall con- 
ditions are imposed at the body. At the outflow, the 
flow is supersonic normal to the boundary and no 
physical boundary conditions are imposed. Further 
details on the boundary treatments can be found 
elsewhere 

The grid used in the Mach = 2.5 case (see Figure 
3), contains 36 radial and 24 circumferential collo- 
cation points. This grid density is sufficient to pre- 
dict a stagnation pressure on the centerline exact 


0-0 n 


9max 9 n 


( 2 ) 


to twelve significant digits. The solutions accuracy 
throughout the domain is monitored using the spec- 
tral decay of polynomial coefficients. The worst por- 
tion of the flow is resolved to eight orders of mag- 
nitude. This solution (referred to from now on as 
the “exact” solution) is spectrally interpolated to a 
sequence of uniformly spaced grids for use as error 
indicators in the finite difference algorithm. Figures 
4, 5 and 6 show the pressure, density and veloc- 
ity vector profiles of the blunt-body flow. Note the 
smooth, well resolved contour lines. 



Figure 5. Density profiles of the “exact” solution 
obtained with the Chebyshev shock fitting algorithm. 


Fourth- Order Linear Scheme 

Conventional high-order central difference dis- 
cretizations lack the robustness necessary to effi- 
ciently converge this problem to steady-state. A 
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Figure 6. Velocity vectors from the “exact” solution 
obtained with the Chebyshev shock fitting algorithm . 

fourth-order upwind biased algorithm, however, 
proves robust and efficient over the Mach number 
range consider in this study. The upwinding uses a 
Lax-Friedrichs splitting of the flux vector. Specif- 
ically, the derivative of the contravariant flux F is 
calculated as: 

DF = \{D+ F+ + D~ F~) 

where F ± = F±X masc U y and D+ = (3,3, 3-4-3), 
D~ = (3 — 4 — 3,3,3). (The nomenclature for deriva- 
tive closure D+ denotes that the left three points are 
treated third-order, the interior is fourth-order and 
the right point is third-order. Details of the nomen- 
clature can be found elsewhere 5 .) The fourth-order 
D + operator is the upwind biased stencil: 

*+i 

U(xi) = ^ aijU(xj) +0(Sx) 4 

j—i- 3 

where row i of the matrix aij is defined as aij = 
-A— [— 1 } 6 , — 18 , 10 , 3 ]. The third-order stencils 
used at all boundary points are the optimal stencils 
derived from nearest neighbor information, biased 
where possible in an upwind direction. The scaling 
parameter A max is the maximum over the entire do- 
main of the contravariant eigenvalue u + c. A similar 
expression is used for the G vector. 

Historically, finite-difference formulations that lo- 
cate a discrete point at the stagnation point are 
susceptible to numerical instability. Careful imple- 
mentation of the boundary conditions eliminates this 
problem. The present formulation is as robust as the 
staggered formulations typically adopted by finite- 
difference algorithms near stagnation points. On 


the inflow plane all variables are prescribed, con- 
sistent with supersonic inflow. On the outflow, the 
solution is simply calculated using one-sided high- 
order information, and no boundary conditions are 
imposed. On the symmetry plane, the v component 
of the velocity is set to zero. On the inviscid wall, 
the no penetration boundary condition is treated 
weakly through the flux. Specifically, the contravari- 
ant flux on the wall is calculated consistent with the 
no-penetration condition. (The normal velocity is 
flipped in sign, to create a reflecting state. The two 
states are then rectified with a Riemann solver to 
obtain a state with no normal velocity). The physi- 
cal velocity at the wall is not required to satisfy the 
no penetration condition. Thus, at steady state, the 
normal velocity at the wall is nonzero, but converges 
to zero with an order property consistent with the 
overall formulation. 

A three stage Runge-Kutta is used to drive the 
solution to steady-state. The calculations are sus- 
pended when the solution error based on the exact 
solution converges to four significant places. 

Third- Order Nonlinear Scheme 

The third numerical scheme employed in this work 
is a third-order accurate ENO algorithm. Spatial ac- 
curacy is achieved by solving the two-dimensional, 
time-dependent Euler equations in control-volume 
form. Complete details of the algorithm are pre- 
sented elsewhere 7 . The adaptive stencils employed 
in the high-order spatial operator are biased in 
smooth regions toward those that are linearly stable 
(See, e.g. 8 ). All required mesh quantities (cell ar- 
eas, grid metrics, etc.) are calculated to consistent 
accuracy 7 . The solution is advanced in time, to- 
wards a steady state, via a three-stage Runge-Kutta 
scheme 9 . Implicit residual smoothing is employed 
in order to accelerate the time integration. 

Fluxes at cell interfaces are approximated in two 
distinct ways. The most accurate treatment is by 
Roe’s approximate Riemann solver, as presented 
in Reference 10. Full machine-zero convergence to 
steady state is not possible with this formulation. 
High-order accuracy is nonetheless attained on the 
grid aligned test case. The robustness and gener- 
ality of the formulation is somewhat questionable, 
making it necessary to use a second formulation. To 
drive residuals to machine zero, it is necessary to 
augment the Roe solver with an entropy fix 11 , and 
to more heavily bias the stencils. These modifica- 
tions give rise to an oscillation near the shock, and 
cause the solution error to revert to first-order ac- 
curacy. We use the straight Roe solver as a means 
of verifying the accuracy of the algorithm, and the 
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second formulation as the general solver used in the 
comparative studies. 

Grids 

The grids used for all calculations originate from 
the Chebyshev spectral solution. A uniform dis- 
cretization in both the radial and circumferential di- 
rections is used. To accomplish shock capturing, the 
grid is extended in the radial direction until exactly 
one quarter of the points are in the supersonic por- 
tion of the flow. The following sequence of grids 
is used to perform a grid refinement study: 33x33, 
65x65, 129x129, 257x257, 513x513. The exact so- 
lution is projected onto these grids (using spectral 
interpolation) for use as the initial condition, and to 
calculate solution accuracies. 


Results 


First-Order motivation 

We begin by presenting a heuristic argument 
explaining why the 2-D steady-state Euler equa- 
tions admit first-order solutions downstream of two- 
dimensional discontinuities. We saw in previous 
works 2 , 3 that time dependent information passing 
though a discontinuity is corrupted to first-order by 
a discontinuity. This statement was quantified us- 
ing the model problem u t + [a(x)u] x = 0 with a(x) 
discontinuous. It was concluded that first-order cor- 
ruption occurs in hyperbolic equations when a(x) 
is discontinuous and information passes through the 
discontinuity. 

We now show that at least part of the steady- 
state 2-D Euler equations are mathematically equiv- 
alent to the model hyperbolic problem, with the 
wave speed a (a) being discontinuous at the captured 
shock. (Here the second spatial dimension takes the 
place of time, resulting in a hyperbolic equation in 
[x,y]). Therefore, based on the 1-D time-dependent 
results, the discontinuous 2-D Euler equations will 
exhibit first-order error downstream of an arbitrary 
captured discontinuity. 

The equations governing the steady-state Euler 
equations in 2-D are 


&F dG 
dx dy 


A dV <9U 
A <9x + B dy 


( 3 ) 


where A = and B = The mathematical 
charactor of equation (3) is governed by the eigenval- 
ues of the matrix M = A -1 B. The eigenvalues of 
the matrix M are A 1j2 = £, A 3)4 = 

For supersonic flow the eigenvalues are real and dis- 
tinct, and the equation set is strictly hyperbolic. For 
subsonic flow, two of the eigenvalues are real and two 


are imaginary and the equation set is mixed elliptic- 
hyperbolic. The two eigenvalues that are real inde- 
pendent of the Mach number are Ai ) 2 = and 
correspond to the advection of information along 
the streamlines -. The mathematical character of 
the hyperbolic portion of the equation is identi- 
cal to that exhibited in the unsteady 1-D equation 
u t + [a(x)u] x = 0 with a(x) discontinuous. Thus, 
the steady 2-D Euler equations will be first-order 
downstream of any captured discontinuity. 

We digress here and note that the multidimen- 
sional accuracy at steady state is different from that 
found in the 1-D case. The 1-D steady state solu- 
tion is design order away from any discontinuities, 
while the 2-D case is first-order behind general dis- 
continuities. We provide no mathematical explana- 
tion for this anomaly other than to say that two 
“dimensions” are necessary to exhibit the first-order 
behavior. 

The Smooth “Near-Body” Problem 

A model problem is needed to test the accuracy of 
the shock capturing algorithms. Different formula- 
tions of the exact solution are used to test the linear 
and non-linear algorithms. We begin by describing 
the testing of the linear algorithm. 

By solving for the flow between the bow-shock 
and the cylinder (the near-body), we test the center- 
line, wall and outflow boundary conditions, as well 
as the design accuracy of the method. Since there 
are no discontinuities in this subproblem, the linear 
formulation should recover design accuracy. Elimi- 
nating the gridpoints outside the bow-shock, a sub- 
sequence of grids is formed having 25x33, 49x65, 
97x129, 193x257, and 385x513 points respectively. 
The boundary conditions for the centerline, wall and 
outflow planes are identical to those described pre- 
viously. The subsonic inflow boundary condition at 
the bow-shock is implemented by solving the Rie- 
mann problem between the numerical state at the 
point and the exact post-shock conditions from the 
spectral algorithm. 

Table la. shows the results of a grid refinement 
study using the linear fourth-order scheme for a 
Mach number of M = 2.5. Reported are the 

I >2 error in density, and the L 2 error in the nor- 
mal velocity component at the inviscid wall, both 
as a fun ction of grid. The £ 2 error is defined as 

L 2 '= y — — ^ , where ^ and are the 
numerical and exact values of the respective func- 
tions at point z, and N is the number of points in 
the norm. The norm of density is formed over the 
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entire domain. (The norm of wall-normal velocity is 
formed over all body points). 



Density 

4 th 

Wall BC 

4 */i 

Grid 

log 1,2 

Rate 

log Li 

Rate 

25x33 

-3.275 


-3.119 


49x65 

-4.230 

3.17 

-4.003 

2.94 

97x129 

-5.256 

3.41 

-4.980 

3.25 

193x257 

-6.330 

3.57 

-6.041 

3.52 

385x513 

-7.419 

3.62 

-7.163 

3.73 


Table la: Smooth test case of near wall region show- 
ing design accuracy of the linear shock capturing 
scheme . 

The convergence rate is asymptoting to the the- 
oretical value of 4. (The third-order closure near 
the boundaries dominates the solution accuracy on 
coarse grids). The normal component of the velocity 
at the wall is going to zero in accordance with the 
weak imposition of the tangency condition. Simi- 
lar convergence behavior is observed for the pressure 
and for other other norms. 

Several specific features of the blunt body, and the 
ENO formulation, allows for a convincing test of the 
non-linear formulation. The full grid is constructed 
with the shock coinciding exactly with an interface 
in the finite volume formulation. If smooth inter- 
polants are used on either side of the cell interface, 
and a Riemann solver which satisfies the Rankine- 
Hugoniot shock jump relations is used, then for this 
special case, the shock is “fit” at the cell interface. 
All these constraints are met by the ENO formula- 
tion using a Roe flux at the interfaces. As noted 
earlier, the temporal residual for this test problem 
would not converge more than four orders of magni- 
tude. In spite of this, design accuracy was achieved 
on the original sequence of four grid. 

Table lb. shows the results of a grid refinement 
study using the third-order ENO scheme for a Mach 
number of M = 2.5. Reported is the L 2 error in 
density, using the same definitions of error reported 
in the linear case. 



Density 

3rd 

Grid 

log £2 

Rate 

33x33 

-3.162 


65x65 

-4.042 

2.92 

129x129 

-5.007 

3.21 

257x257 

-6.045 

3.45 


Table lb: Special ^fitted” test case showing design 
accuracy of ENO formulation using the straight Roe 
flux at interfaces . 

The convergence rate is asymptoting to a value 
greater than the theoretical accuracy of 3. Note in 



Figure 7. Density profiles of the captured solution 
obtained with the third-order ENO method \ using 
only the Roe flux at the interfaces. The solution is 
third- order accurate. 


an absolute sense, that the errors from the third- 
order ENO formulation are comparable with those 
obtained from the fourth-order linear algorithm. 
Two dimensional plots of the pressure and density 
obtained from the ENO algorithm are nearly indis- 
tinguishable from the “exact” solutions shown in 
Figures 4 and 5. Figure 7 shows the pressure profiles 
obtained with the nonlinear ENO algorithm for this 
test case. Note the monotone one point bow-shock. 

Several points can be deduced from these studies. 
First, the shock capturing algorithms exhibit design 
accuracy on the near-body subproblem. This indi- 
cates that the symmetry, wall and outflow bound- 
aries are correct. Secondly, the agreement to seven 
significant digits between the capturing algorithm 
and the Chebyshev spectral solution indicates that 
the “exact” answer is correct to at least 7 digits. 

The Shock-Captured “Blunt- Body” Problem 

We now use the validated linear fourth-order and 
the nonlinear third-order algorithms to capture the 
bow-shock around a Mach 2.5 circular cylinder. The 
only modifications necessary to capture the blunt- 
body shock are, 1) in the linear formulation, the 
original grids are used, and supersonic free-stream 
conditions are imposed at the inflow. 2) For the 
nonlinear ENO formulation the Roe solver at the 
interfaces is modified with an entropy fix, and the 
stencil biasing parameters were tuned. These mod- 
ifications allow a steady-state residual to be driven 


6 






to machine precision, but eliminate the possibility of 
“fitting” the shock at the interface. 
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Figure 8. Finite difference grid employing 65x65 
evenly spaced points, used to capture the Mach = 
2.5 blunt-body flow . 
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Figure 9. Pressure profiles of the captured solution 
obtained with the fourth- order finite difference algo- 
rithm. 

Figure 8 shows the 65x65 grid used in both formu- 
lations. Figures 9 and 10 show the two-dimensional 
pressure and density profiles obtained with the linear 
method. Figure 11 shows the pressure profiles ob- 
tained with the nonlinear method. All are obtained 
on the 65x65 grid. The solution from the linear al- 
gorithm has huge oscillations at the captured shock, 
whereas the nonlinear formulation is nearly mono- 
tone. 

Refinement studies using the original sequence of 
grids are presented in Tables Ha and lib. In table 



Figure 10. Density profiles of the captured solution 
obtained with the fourth-order finite difference algo- 
rithm. 

Ha, the L*i norm of the density and pressure errors 
are presented, as obtained with the fourth-order lin- 
ear algorithm. The density and pressure norms are 
formed using only the half of the domain closest to 
the body, in order to exclude any points close to the 
captured shock region. 



Density 

4^ 

Pressure 

4 th 

Grid 

log £2 

Rate 

log L>2 

Rate 

33x33 

-1.756 


-1.651 


65x65 

-3.404 

5.47 

-3.370 

5.74 

129x129 

-3.708 

1.01 

-3.778 

1.32 

257x257 

-4.014 

1.02 

-4.117 

1.13 

513x513 

-4.317 

1.01 

-4.437 

1.06 


Table Ha: Captured Mach 2.5 bow-shock around cir- 
cular cylinder using linear shock capturing scheme. 

In marked contrast to the accuracy of the near- 
body problem, the 2-D captured shock accuracy is 
asymptotically first — order in space. On coarse 
grids the scheme converges near the design accuracy, 
but is dominated by the first-order shock error on 
finer grids. Note that the density and pressure ex- 
hibit nearly identical behavior. We note, but do not 
show that the wall penetration error is asymptoting 
to design accuracy, as it did in the near-wall test 
problem. 

In table lib, the the L 2 norm of the density er- 
rors are presented, as obtained with the third-order 
ENO algorithm. The same error norms are used as 
in the linear algorithm. The 2-D captured shock ac- 
curacy is asymptotically first — order in space. The 
nonlinear formulation yields increased accuracy on 
coarse grids, but surprisingly is less accurate asymp- 
totically than the fourth-order linear algorithm. Not 
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Figure 11. Density profiles of the captured solution 
obtained with the third- order ENO algorithm , using 
the Roe flux with an entropy fix at the interfaces . 
The solution is first-order accurate . 

surprisingly, there is no component of high-order 
convergence on this problem, as there is with the 
linear scheme. Work continues with the nonlinear 
scheme to try to recover better asymptotic error be- 
havior at steady state. Based on previous experience 
with the nonlinear formulation, the asymptotic error 
level can be reduced by further optimization of the 
biasing parameters and the entropy fix. 



Density 3 rd ENO 

Grid 

log Li Rate 

33x33 

65x65 

129x129 

257x257 

-2.275 

-2.572 0.99 

-2.871 0.99 

-3.170 0.99 


Table lib: Captured Mach 2.5 bow- shock around cir- 
cular cylinder using third-order ENO shock captur- 
ing scheme . 

Comparing the “near-body” and “blunt-body” so- 
lution errors, we have established that capturing the 
discontinuity in the blunt-body case causes a degra- 
dation in accuracy. The resulting convergence rate 
is first-order for both the linear and nonlinear for- 
mulations. We claim that this result is general for 
any multi-dimensional shock, and any high-order nu- 
merical method. A heuristic explanation for this 
phenomena is related to the ambiguous shock po- 
sition within the cell. Capturing the shock provides 
(through conservation and the Rankine-Hugoniot re- 
lation) the post shock conditions for the smooth 


near-body problem. Because the exact position of 
the shock is ambiguous to order Sx these “numer- 
ical” conditions are incorrect to first-order. This 
first-order error is relatively small on coarse grids, 
but at some level of refinement becomes the domi- 
nant error term. 

Possible exceptions to the first-order conclusion 
include the case where a multi-dimensional shock 
can be “one-dimensionalized” by a suitable rotation. 
For example, if the numerical grid is aligned with the 
shock throughout the domain, and the Roe solver is 
used to rectify the interface states then design ac- 
curacy is possible at steady-state. This is roughly 
equivalent to shock fitting, and is a formidable task 
in general situations. 

First-Order Implications 

We have demonstrated that captured shocks in 
2-D are first-order, regardless of the design accu- 
racy of the capturing scheme. At some point of 
spatial resolution the solution is dominated by the 
first-order component of the error. This raises the 
obvious question, “If solution error is bounded at 
first-order, to what extent should high-order for- 
mulations be considered for capturing discontinu- 
ous flows”? Specifically, for what class of problems 
does the first-order error component dominate the 
solution error, making additional high-order work 
counter-productive? We first present a parametric 
study of design accuracy, to quantify the nature of 
the first-order error component. We then begin to 
address the specific role of high-order algorithms in 
the context of shock capturing. 

Dependance of Error on Design Accuracy 

Tables Illa-IIIc contain a comparison of differ- 
ent design accuracies on the Mach 2.5 blunt body 
problem. The study is performed using the linear 
algorithm. Work continues with the nonlinear al- 
gorithm to quantify the generality of the conclu- 
sions drawn in this study. We include spatial op- 
erators of first- second- and fourth-order accuracy. 
The fourth-order pressure and density results shown 
previously are compared with those obtained with 
the second- and first-order algorithms. We see that 
all methods are asymptotically first-order. The re- 
sults from the first-order (design accuracy) scheme 
are very poor, and only a partial listing is included. 
Note that the asymptotic level of error in the second- 
and fourth-order schemes is the same. It appears 
that the first-order component of the captured er- 
ror is almost independent of the design order of 
accuracy of the method. The fourth-order scheme 
quickly approaches the asymptotic limit, while the 
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second-order scheme asymptotes more slowly. The 
first-order (design accuracy) scheme will never be 
dominated by shock error. 



Density 

4 th 

Pressure 

qt/i 

Grid 

log £2 

Rate 

log L 2 ■ 

Rate 

33x33 

-1.756 


-1.651 


65x65 

-3.404 

5.47 

-3.370 

5.74 

129x129 

-3.708 

1.01 

-3.778 

1.32 

257x257 

-4.014 

1.02 

-4.117 

1.13 

513x513 

-4.317 

1.01 

-4.437 

1.06 


Table Ilia: Fourth-order linear algorithm . 



Density 

2 nd 

Pressure 

2 nd 

Grid 

log 1/2 

Rate 

logI 2 

Rate 

33x33 

-2.327 


-2.087 


65x65 

-2.776 

1.49 

-2.933 

2.81 

129x129 

-3.242 

1.55 

-3.664 

2.43 

257x257 

-3.660 

1.39 

-4.321 

2.18 

513x513 

-4.034 

1.24 

-4.452 

0.44 


Table Illb: Second- order linear algorithm . 



Density 


Pressure 

l$t 

Grid 

log £2 

Rate 

log 1-2 

Rate 

33x33 

-0.891 


-0.745 


65x65 

-1.335 

1.47 

-1.049 

1.01 

129x129 

-1.683 

1.16 

-1.348 

0.99 

257x257 

-2.040 

1.19 

-1.652 

1.01 


Tables IIIc: First-order linear algorithm used to cap- 
ture the Mach 2.5 bow-shock around a supersonic 
circular cylinder. 

A heuristic model expressing the nature of the so- 
lution error exhibited by the linear schemes, must 
include a first-order and a design-order component. 
A model approximately consistent with the solution 
errors presented in Tables Illa-c is: 


£ = Cl ((fa) 1 + c 2 (<fa) r (4) 

where e is the solution error, r is the design order of 
the numerical algorithm, and c\ and C2 are problem 
dependent constants. For r > 1, and any finite val- 
ues ci and C2, the solution error will asymptotically 
be dominated by the first-order error component. If 
ci << C2j the solution will exhibit high-order con- 
vergence on coarse grids. 

The Role of High- Order Shock Capturing 
Practitioners, are interested in absolute error, not 
the order of the spatial approximation. Order of ac- 
curacy and absolute error are closely related in the 


CONVERGENCE : MACH = 2 



Log (dx) 

Figure 12. Grid refinement study at Mach = 2.0 , 
showing dependence of solution error on the design 
accuracy of the numerical scheme . 

asymptotic limit of high accuracy. Most engineer- 
ing problems, however, require error levels far re- 
moved from this asymptotic limit. Thus, a second- 
order algorithm could converge at its design rate, 
for many target accuracies. This would explain why 
no one has quantified the first-order convergence of 
captured shocks over the past thirty years. 

Practitioners are beginning to demand more ac- 
curacy from their calculations. It is now common 
in aero-acoustics and electro-magnetics to demand 
resolution of scales which differ by several orders of 
magnitude. In configuration aerodynamics a com- 
monly used metric of solution accuracy for viscous 
airfoils is one drag count, which corresponds to so- 
lution accuracies of at least four significant digits. 
As the design accuracy requirements increase, it is 
much more likely that the solution error will be dom- 
inated by the first-order component downstream of 
a discontinuity. When this occurs, global refinement 
of the solution is no longer a viable procedure; local 
refinement of the shock region must be performed to 
obtain an accurate solution. 

We now show a parametric study in Mach num- 
ber, using the linear algorithms of fourth-, second- 
and first-order accuracy. The objective of the study 
is to show that the efficacy of high-order methods 
is dependent on the problem as well as the target 
accuracy level. Figure 13 shows a plot of the pres- 
sure data presented in Table Ilia. Lines of unit slope 
are included to compare the asymptotic convergence 
rates of the methods. On this plot one can more eas- 
ily identify which method is the most efficient for a 
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CONVERGENCE : MACH = 2.5 



Log (dx) 

Figure 13. Grid refinement study at Mach = 2.5, 
showing dependence of solution error on the design 
accuracy of the numerical scheme . 

given target accuracy. At a target accuracy of 10“ 2 
the second order scheme is most desirable. (On the 
coarsest grid the second-order scheme is more ac- 
curate than the fourth-order scheme.) At a target 
of 10“ 3 the fourth-order scheme is most desirable. 
At one drag count, or four significant digits (10“ 4 ), 
both methods are in the first-order asymptotic limit. 
Neither method is computationally efficient, and re- 
finement of the grid in the bow-shock region would 
be the best use of computer resources. 

The results in figure 13 are consistent with the ob- 
servation that in many circumstances the first-order 
shock component of the solution error is negligibly 
small. Note that the second-order scheme is pro- 
foundly more accurate than the first-order scheme, 
and only begins showing its first-order character on 
extremely fine grids. Seldom have accuracy levels 
comparable to those reported here been obtained by 
CFD practitioners because of the costs involved. It 
is not surprising that the first-order nature of the 
captured solution has not been uniformly recognized 
till now. 


Mach 

Density 

Pressure 

2.0 

2.667:1 

4.500:1 

2.5 

3.333:1 

7.125:1 

3.0 

3.857:1 

10.333:1 


Table IV: Centerline, bow-shock, density and pres- 
sure ratios, as a function Mach number . 

Figures 12 and 14 show results for Mach numbers 
2 and 3, respectively. The three numerical schemes 
used in the Mach = 2.5 case are used in these stud- 
ies. Table IV shows the shock strengths as a function 


CONVERGENCE : MACH = 3 



Log (dx) 

Figure 14. Grid refinement study at Mach = 3.0, 
showing dependence of solution error on the design 
accuracy of the numerical scheme. 

of Mach number for the range 2 < M < 3. The 
conclusions drawn in the Mach = 2.5 case appear to 
be reasonably general for this range of Mach num- 
bers. Specifically, the first-order (design accuracy) 
results are not competitive with either the second- 
or fourth-order results. Both higher-order methods 
exhibit design accuracy on coarse grids, but asymp- 
tote to first-order on fine grids. The asymptotic level 
of the higher-order schemes is nearly independent of 
the design accuracy. At low target accuracies, the 
second-order scheme is most desirable, and at mod- 
erate target accuracies, the fourth-order scheme is 
most desirable. At high target accuracies, the high- 
order schemes are first-order because of the shock, 
and regridding the shock is necessary to effectively 
decrease the overall accuracy of the calculation. 

The convergence rates at all Mach numbers are 
consistent with the heuristic model proposed in 
equation 4. In addition, we see a slight depen- 
dence of the asymptotic first-order shock error on 
the shock strength, that increases with Mach num- 
ber. Work continues to quantify this effect, and to 
be able to predict the shock capturing error levels at 
the shock. It is hoped that this quantification will 
provide some general guidance for capturing shocks 
of different strengths, and for when refinement of the 
shock is necessary to ensure satisfactory results near 
the body. 

Conclusions 

An assessment of the accuracy of shock capturing 
schemes is made for two-dimensional steady flows. 





A linear fourth-order method and a nonlinear third- 
order method are used in this study. It is shown, 
contrary to conventional wisdom, that captured two- 
dimensional shocks are asymptotically first-order, 
regardless of the design accuracy of the numerical 
method used to obtain them. 

The solution error of the captured discontinu- 
ity can be characterized by a first-order component 
ci(Sx) 1 , and a high order component C 2 (Sx) r (r is 
the design order of the method). With sufficient 
resolution, the error is dominated by the first-order 
terms. In a practical sense, however, frequently the 
solution error is dominated by the high-order com- 
ponent, giving the appearance of a high-order con- 
vergence rate. 

Finally, we acknowledge that demonstrating a 
first-order convergence rate with two numerical 
methods, does not prove that captured shocks are 
destined to be first-order accurate. Rather, it 
demonstrates that we do not presently know how 
to achieve high order convergence in a general con- 
text when discontinuities are present. Work contin- 
ues towards the goal of design accuracy for captured 
discontinuities. 
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ABSTRACT 

The numerical study of aeroacoustic problems places stringent demands on the choice of a 
computational algorithm, because it requires the ability to propagate disturbances of small amplitude 
and short wavelength. The demands are particularly high when shock waves are involved, because 
the chosen algorithm must also resolve discontinuities in the solution. The extent to which a 
high-order-accurate shock-capturing method can be relied upon for aeroacoustics applications that 
involve the interaction of shocks with other waves has not been previously quantified. Such a 
study is initiated in this work. A fourth-order-accurate essentially nonoscillatory (ENO) method 
is used to investigate the solutions of inviscid, compressible flows with shocks. The design order 
of accuracy is achieved in the smooth regions of a steady-state, quasi-one-dimensional test case. 
However, in an unsteady test case, only first-order results are obtained downstream of a sound- 
shock interaction. The difficulty in obtaining a globally high-order-accurate solution in such a case 
with a shock-capturing method is demonstrated through the study of a simplified, linear model 
problem. Some of the difficult issues and ramifications for aeroacoustic simulations of flows with 
shocks that are raised by these results are discussed. 
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INTRODUCTION 


This work is motivated by the desire to develop numerical methods that will be useful in the 
study of aeroacoustic phenomena that occur in flows with shocks. For example, shocks in jet 
flows, on wings, and in supersonic combustion inlets contribute significantly to sound generation. 
Problems such as these represent some of the more challenging aspects of ongoing research in the 
developing area of computational aeroacoustics (CAA). 

One of the purposes of this work is to open a discussion on the relative merits of the numerical 
methods which can simulate sound sources that are generated in flows with shocks. For a compu- 
tational algorithm, obtaining acoustic information from a numerical solution that involves shock 
waves is a demanding proposition. In general, high-order accuracy is required for the propagation 
of high-frequency low-amplitude waves. In addition, the shock must be adequately captured. A 
significant body of work has been devoted to the development of numerical schemes that possess 
both properties. 

The contribution of a shock wave to the generation of sound is attributed to the motion of the 
shock wave and its interaction with other disturbances in the flow. These disturbances can be large, 
such as changes in the mean flow, or small, such as acoustic or vortical disturbances. The extent to 
which such interactions must be accurately predicted in order to propagate reliable information to 
an observer will be investigated. 

In the following section, some of the currently available methods that vary in their approach to 
high-order accuracy and shock capturing are surveyed. Next, a fourth-order-accurate essentially 
nonoscillatory (ENO) method is applied to the steady solution of a nozzle flow with a shock, 
and fourth-order-accurate results are obtained. However, in the application of the same numerical 
scheme to a time-dependent problem, that of a sound-shock interaction, accuracy suffers. The 
disappointing results in regard to the accuracy of this solution are explained through the study 
of a simpler linear model problem. The numerical scheme is modified with a subcell resolution 
technique in order to obtain a globally high-order-accurate solution for time-dependent simulations. 
The original sound-shock interaction problem is again solved with the modified scheme, and 
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a fourth-order-accurate solution is achieved. Some of the difficult issues and ramifications for 
current methods that are raised by these results are discussed in the final section. 

HIGH-ORDER ACCURATE SHOCK-CAPTURING METHODS 

The high-order-accurate spatial operator that is desired in a shock-capturing method for CAA 
problems requires that information be taken from a large number of discrete locations within the 
solution. Such an operator will cause large oscillations in a discontinuous solution unless special 
precautions are taken. Many methods are available in the literature that attempt to balance the 
properties of high-order accuracy and shock capturing. However, they can be classified into two 
categories i.e., linear and nonlinear. 

Within the linear class of numerical shock-capturing schemes, the interpolation set for the 
approximation of the solution or its derivatives is fixed as a function of grid location. Linear methods 
admit oscillations in regions in which physical gradients are inadequately resolved. Central- 
differencing operators and spectral methods are particularly prone to these numerical oscillations. 
For nonlinear problems, limiters or filters are usually required to keep oscillations from growing 
without bound. The general features of these ideas are demonstrated in the work of Don , 1 in which 
pseudospectral methods are used for compressible flow problems with shocks by locally applying 
a simple filter to keep the solution bounded during the computation. The final solution is post- 
processed in order to remove the oscillatory information that developed in time. Nothing in this 
approach mandates the use of pseudospectral techniques; other higher-order schemes (e.g., central- 
difference or compact schemes 2 ) could also be used. For example, Harten and Chakravarthy 3 
suggest a polynomial interpolation procedure in which the coefficients of the polynomial are 
limited by a switch that makes a coefficient very small if the corresponding solution derivative is 
discontinuous. 

In the nonlinear class of schemes, the strategy with respect to discontinuities is to employ some 
sort of adaptive interpolation. The goal is to achieve formal high-order accuracy in smooth regions 
and high shock resolution without oscillations. The class of ENO schemes 4,5 has been designed to 
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have such properties. As originally presented, the local polynomial approximation operator adapts 
its interpolation set to the smoothest available part of the solution. This adaptive interpolation 
strategy has also been employed in the development of psuedospectral hybrid schemes. 6-8 

Although discontinuous solutions generated by a linear strategy are usually not as pictorially 
pleasing as solutions in which shock profiles are monotone, linear schemes are more computa- 
tionally efficient than nonlinear schemes. The efficiency of a numerical algorithm is extremely 
important for aeroacoustic simulations because such problems are time dependent and require a fine 
computational mesh for the resolution of high-frequency disturbances. Because nonlinear methods 
are designed to avoid the production of spurious oscillations, the stability of a calculation of a flow 
with shocks is more readily obtained. However, their adaptive interpolation operator significantly 
hampers their efficiency relative to linear schemes. 

In the literature to date, the design accuracy of a numerical method, whether linear or nonlinear, 
has been demonstrated with solutions to smooth problems. Problems with discontinuous solutions 
are most commonly used only to illustrate the ability of the scheme to obtain a stable solution to 
such problems. In the following two sections, the ability of a fourth-order accurate ENO method 
to achieve high-order accuracy in the smooth regions of a flow with a shock is investigated. The 
method is applied to obtain the solutions of a steady problem and a time-dependent sound-shock 
interaction. The accuracies of these solutions are then analyzed through the study of a simpler 
linear model problem. 


STEADY SHOCK IN A NOZZLE 

A steady-state flow with a shock in a quasi-one-dimensional converging-diverging nozzle is nu- 
merically investigated. The governing equations are the quasL-one-dimensional Euler equations: 

+ r^ AF) = G (la) 


where 



' p ' 


pu 


0 1 

u= 

pu 

pE 

, F = 
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The variables p,u,P,E , and A are the density, velocity, pressure, total specific energy, and nozzle 
area, respectively. The equation of state is 


P = ( 7 -l )p{E - ^ u 2 ) 


where 7 is the ratio of specific heats, which is assumed to have a constant value of 1.4. 

The spatial domain of the nozzle is 0 < x < 1 , and the flow is oriented from left to right. 
The nozzle shape is determined exactly through the requirement of a linear distribution of Mach 
number from M = 0.8 at the inlet to M — 1.8 at the exit, assuming the flow is isentropic and fully 
expanded. The resulting area distribution A(x) is illustrated in Fig. la. The flow variables are 
normalized with respect to upstream stagnation conditions and the area with respect to the value at 
the throat, x = 0.2. 

Given the prescribed area distribution, the Mach 0.8 inflow state is retained at x = 0, and the 
outflow condition at x = 1 is determined such that a shock forms at x s = 0.5, which corresponds 
to a preshock Mach number of M — 1.3 . A steady-state solution is obtained by implementing a 
fourth-order-accurate ENO method until residuals are driven to mac hi ne zero. Spatial accuracy is 
achieved by solving the equations in control- volume form as presented in Ref. 4. The equations are 
integrated in time via a third-order-accurate Runge-Kutta scheme. 5 This numerical method will be 
referred to as “ENO-4-3.” As has been established in previous research, 9-12 the adaptive stencils 
employed in the spatial operator are biased in smooth regions toward those that are linearly stable. 
Fig. lb depicts the solution, with respect to Mach number, on a uniform mesh of 128 cells. 

One of the simpler methods of determining the error of this solution relies on the fact that the 
value of the entropy-like quantity Pj p 1 is piecewise constant: 


S = 


P 

P 1 



X < x s 
X > x s 


The subscripts of S denote the pre-shock and post-shock stagnation values, respectively. This 
quantity is plotted in Fig. lc. The pointwise entropy error for this solution on four successively 
refined meshes is illustrated in Fig. 2. Clearly, the accuracy is fourth order on either side of the 
shock, as demonstrated by the error data in Table I. The variable N c is the number of cells. The 
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errors ||e|| are computed in the Li and norms. The number r c is the computational order of 
accuracy and is determined by the slope of these tabulated values on a log-log plot: 

c In (hk/ht+x) 

where is the error measured on a mesh of uniform spacing hk with Ti*, > hk+\ for k— 1,2, 3. 

Although these results are encouraging, the time independence of the solution makes this a 
convenient example for the demonstration of high-order accuracy in the presence of a shock. As 
will be shown in the study of an unsteady problem in the following section, a moving shock presents 
a greater challenge in regard to high-order-accurate shock capturing. 

SOUND-SHOCK INTERACTION 

The interaction of a sound wave with a shock in a one-dimensional flow is numerically in- 
vestigated. The effects of shocks on sound waves, and vice versa , are important to the acoustics 
and performance of aircraft design. Therefore, the ability to obtain an accurate solution to such a 
model initial-boundary-value problem (IB VP) is important in the development of shock-capturing 
methods for CAA research. Similar one-dimensional problems have been the subject of other 
studies. 13,14 

The governing equations are the one-dimensional Euler equations: 

^ + 1^) = 0 (2a) 

where the components of U and F(U ) are identical to those given in Eq. lb. The equation of state 
is also the same as in the previous example. 

The spatial domain is 0 < x < 1 . The piecewise constant initial conditions, Ul and Ur, are 
those of a steady shock located &tx s = 0.5. The flow is from left to right, and the state Ul is a Mach 
2 flow upstream of the shock. The flow variables are normalized with respect to this upstream flow. 
At t = 0, an acoustic disturbance is introduced at x = 0: 
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(2b) 


P(0,t) 

p{o,t) 

•u(0, t) 


Pl ( 1 + e sin ut ) 
1 h 


= Pl 


= u L + 


Pl 

2 


7 - 1 


c(0,£) - c L 


where u is the circular frequency, e is the amplitude, and c = y^Pj p is the local sound speed. 

The numerical solution of this problem is obtained through the implementation of the ENO-4-3 
algorithm. This algorithm will be fourth-order accurate even for a time-dependent problem when 
the time step is suitably restricted. The exact solution is obtained by a two-domain Chebyshev 
spectral technique. 15 Shock fitting is used to divide the domain into two computational regions. A 
Chebyshev collocation method is used in each region for the spatial discretization. A fourth-order 
Runge-Kutta scheme is used to discretize time. Sufficient spatial and temporal resolution are used 
to guarantee machine precision of the solution. 

Fig. 3 depicts the pressure perturbation 6P(t ) = P(x, t) — P(x, 0) at t = 30 T\, where 
T\ — 27r/u; is one period of the incoming acoustic wave. The acoustic wave amplitude is e = 0.001, 
and u = 2irk(uL + cl) is determined by requiring a wave number Jc = 6 with respect to unit length 
and a mean wave speed ui + cl- The calculation, represented by circles, was performed on a 
uniform mesh of 256 cells with a Courant number of 0.5. The exact solution is represented by a 
continuous line. In this pictorial measure, the numerical algorithm performs well with respect to its 
prediction of the amplified sound wave at higher frequency downstream of the shock. The missing 
circle values near the shock are off the plot and are due to the use of the stencil-biasing parameters 
near a moving discontinuity. 16 

Even more instructive, however, is the pointwise error made by this calculation with respect to 
the mesh width. Fig. 4 illustrates this error on four successively refined meshes. The solution is 
clearly fourth-order accurate upstream of the shock, but only first-order downstream of the shock, 
as shown by the L error data in Table II. The errors are computed on two spatial subdomains: 
0 < x < 0.45 and 0.55 < x < 1. In this manner, the first-order error that is generated in the 
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vicinity of the shock is avoided. This disappointing result is more easily explained through the 
study of a simpler model problem, which is examined in the following section. 


A LINEAR MODEL PROBLEM 


The lower-order accurate results of the previous section can be analyzed through the study of a 
linear scalar model problem. The solution of the following IB VP is instructive because it isolates 
the important phenomenon of propagation of information through a discontinuity. This trait will 
be common to almost any aeroacoustic problem that involves shock waves. Consider the scalar 
equation: 


du , .du 

w + “ (a % = 0 


(3a) 


where the piecewise constant wave speed a(x) is 


a(x) = 


2 , x < x s 
1 , X > x s 


(3b) 


The initial conditions are chosen as 


u(x, 0) 


I , X < x s 
1 , X > x s 


(3c) 


The spatial domain is 0 < x < 1 , and the discontinuity location is x s = 0.5. The inflow boundary 
condition is 

u(0, t) = i ( 1 + e sin u>t ) (3d) 

with e = 0.001 and u determined by requiring a wave number k = 2 with respect to unit length 
and the upstream wave speed a = 2. 

For the purpose of the subsequent discussion, it is necessary to briefly describe the numerical 
scheme. The semi-discrete, finite-volume formulation is obtained by integrating Eq. 3a on an 
interval [x,_ 1 / 2 , x,- +1 / 2 ] with center x, and "volume" A x; : 



-1 ’ 

Ax,- . 




(4a) 
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where 


_ . . 1 r x i+ 1/2 

u i\t) = T — / u[x,t)dx 

is the cell average of u on the zth interval at time f , and the flux fi + i/ 2 (t) is 


(4b) 


fi+i/ 2 (t) = a(x {+1/2 )u(xi +l/ 2 ,t ) 


(4c) 


Temporal integration is achieved by a Runge-Kutta method. 5 A numerical approximation to the 
flux in Eq. 4c is determined in two steps. First, given the cell averages, the solution is approximated 
pointwise within each cell. 


Vi(x) m u(x, t) , Xi-1/2 < x < Xi+ 1/2 (4d) 

In this particular application, Vi( x) is a cubic polynomial. The process by which Vi(x) is obtained 
will be referred to as “reconstruction.” This reconstruction operator contains the adaptive stencil 
algorithm, which avoids interpolation across steep gradients. (See Ref. 4 for details.) Then, at a 
cell interface x t - +1 / 2 , two solution values can exist, as shown in Fig. 5. Correspondingly, two flux 
values can also exist. The second step, then, is simply to choose the upwind value which, for the 
present problem, yields the following numerical flux approximations: 

fi- i/a(<) » , 4e , 

/*+ 1/2M » a(x i+ 1 / 2 )Vi(xi +1/2 ) 

Fig. 6 depicts the perturbation 6u(t) = u(x,t ) — u(x, 0) at t = 10 T\, with the ENO-4- 
3 algorithm. Note the similarity of the features of this solution to those of the sound-shock 
interaction in regard to the changes in amplitude and frequency across the discontinuity. The 
calculation, represented by circles, was performed on a uniform mesh of 64 cells with a Courant 
number of 0.5. The exact solution is 


f | [ 1 + e sin u>(t — x/2) ] , 0 < x < x s 

1 1 + e sin u(t — x + x s ) , x s < x < 1 


( 5 ) 


for t > 3/4, after the perturbation first reaches x = 1. 

Because of the linearity of this problem, the discontinuity location x s remains fixed for all 
time, unlike the sound-shock interaction problem. Therefore, x s can be conveniently placed with 
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respect to a given mesh. For instance, if the computational mesh is composed of an even number of 
uniform cells, then x s = 0.5 will lie on a grid point (i.e., a cell face). Fig. 7a depicts the pointwise 
computational error for this problem on four successively refined meshes for which x s lies on a 
grid point. The numerical solution is fourth-order accurate throughout the domain, as shown by 
the error data in Table III. However, for an odd number of uniform cells, x s is interior to a cell, and 
the solution is first-order accurate downstream of the discontinuity, as shown in Fig. 7b and Table 
IV. 

The strikingly different results in Figs. 7a and 7b yield insight into the disappointing results for 
the sound-shock interaction problem of the previous section. Even if the initial steady shock were 
located on a cell face, it would move into the interior of a cell upon interaction with the upstream 
acoustic wave. Therefore, the following explanation for the lower-order accurate results focuses 
on the location of the discontinuity with respect to the boundary or interior of a cell. 

The difference in the results in Figs. 7a and 7b can be explained in terms of the interpolation 
operator and its direct influence upon the numerical flux. First, consider the case in which the 
discontinuity is on a grid point, for example, the left-hand endpoint of the zth cell, as shown in Fig. 
8a. Because of the adaptive interpolation, the polynomial approximation within the adjoining cells 
is of pointwise, high-order accuracy. The large jump at the left-hand interface of the ith cell is 
immaterial because the flux is determined by the value Pi-x^t-i/a) an d the upstream wave speed 
a = 2, as given by Eq. 4e. Now, consider the case in which x s is in the interior of the ith cell, as 
in Fig. 8b. Unlike Fig. 8a, the polynomial approximation in the ith cell now contains a pointwise 
error that is 0(1). This still does not affect the flux into the ith cell, because the inbound flux is 
determined only by information in the ( i — l)th cell. However, the outbound flux is influenced 
by the value 7 ? ,(x, + 1 / 2 ), which is ultimately responsible for the large errors downstream of the 
discontinuity. 

The loss in accuracy in numerical solutions of linear problems with discontinuous initial data 
has been the subject of previous research by other authors. 17-20 All of these previous studies 
involved solutions of coupled linear systems. It is, therefore, instructive to note that the solution of 
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IB VP 3 is analogous to that of a coupled system in the following way. Let ui denote the solution of 
Eqs. 3a-d on 0 < x < x s , and u 2 denote the solution on x s < x < 1. Now consider, for 6 > 0, the 
solution u 2 on x s < x < x s + 6. The downstream solution u 2 is coupled to u x by the requirement 
of flux conservation: 

lim a(x s + 6) u 2 (x s + 6,t) — a(x s ) Ui(x s , t) 

6 —*- 0 

As in the present study, the work of Majda and Osher 17 was concerned with the inherent degradation 
in accuracy in a region in which information must be numerically propagated across a discontinuity. 
Majda and Osher suggested that this difficulty could be circumvented by smoothing the initial data. 
(See Ref. 17 for details.) However, because the goal of the present work is the application of 
numerical methods to solutions of nonlinear problems, the approach suggested by Donat and 
Osher 20 is more appropriate to use here. These authors propose to maintain accuracy across 
discontinuities by using Harten’s ideas on subcell resolution. 21 

In an attempt to achieve a globally high-order-accurate solution, a subcell resolution technique 
is considered. In particular, a modification of the procedure presented in Ref. 21 is detailed. The 
goal is to obtain a better pointwise approximation of the outbound flux in the ith cell in Fig. 8b. 
If the ith cell is decomposed into two subcells whose adjacent face is x s , then the polynomial 
approximation Vi(x) would be replaced with the piecewise polynomial approximation 


f Vl{x) , Xi_ 1/2 < X < X s 
\ Vr(x) , X s < X < X i+ i / 2 


( 6 ) 


where Vl{x) and Vr(x) approximate u(x, t) in the ith cell to high-order pointwise accuracy on 
either side of the discontinuity. Because the inbound flux is accurately approximated with the 
polynomial in the ( i — l)th cell, Vl(x) can be determined by simple extrapolation of Vi-i(x), as 
shown in Fig. 9. With this reasoning, for the present piecewise cubic polynomial reconstruction, 


clearly 

1 r x s r x i+i/2 

Ui(t) = — — / Vi-i(x) dx + / Vr(x)cIx + 0(Aa; 4 ) (7) 

1/2 ^ x 3 

Finally, Vr(x) must be accurately determined, which enables an accurate outbound flux cal- 
culation given by 'Pr(x; + 1 / 2 ). Let ur denote the subcell average of Vr(x) on [x s , a; i+1 / 2 ]. Then, 
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clearly, from Eq. 7, 

+ 0( Ax 4 ) (8) 

The desired cubic polynomial Vr(x) is then determined by applying the reconstruction operator to 
the set {«/{, iii + i(t), u l+2 (i), ^+3(0}’ which contains only smooth solution information. 

Note that, for the solution of the present problem, the outbound flux cannot be determined by 
the leftward extrapolation of Vi+ i(x), as in Ref. 21. This method would be appropriate only if 
the initial conditions contained the propagating wave profile. In the current problem, at t = 0, the 
solution is constant downstream of x s . Therefore, the extrapolation of V l+ i(x) into the zth cell 
would never allow the wave to propagate into the (i + l)th cell. 

The subcell resolution method described above was incorporated into the ENO algorithm, and 
the linear model problem was again solved with the discontinuity located in the interior of a cell. 
This modified scheme will be referred to as “ENO-4-3-SR.” The globally high-order-accurate 
results are shown by the pointwise error plot in Fig. 10 and the data in Table V. 


Ur = 


x 


t+l/2 


— X 


J rx 3 

' Vi-\(x) dx 

x i- 1/2 


SOUND-SHOCK INTERACTION: REVISITED 


The sound-shock interaction problem is again solved; this t im e the subcell resolution methodol- 
ogy from the previous section is incorporated. However, the subcell technique applied in the linear 
problem above does not immediately carry over to this nonlinear case. The method in Eqs. 6-8 is 
dependent upon the discontinuity location x s . For the nonlinear case, at a given time t, the shock 
location is not known a priori and, therefore, must be determined before the subcell technique in 
Eqs. 6-8 can be applied. This step can be accomplished by solving H(x s ) = 0 for x s in the interval 
(xi- 1 / 2 , ^i+ 1 / 2 ) where the function H{x s ) is given by 


H(x.) = 


A Xi 


r x > r x i+ 1/2 

/ Vi~i(x)dx + / Vi+i(x) dx 
J x s 


- Ui(t ) 


( 9 ) 


Clearly, this function is derived from the relationship in Eq. 7. 

For the current ENO-4-3-SR algorithm, H(x s ) is a quartic polynomial. For a sufficiently small 
A Xi, a unique real root of H(x s ) exists in the interval (a; J _ 1 / 2 , X{+i/ 2 ). The root of H(x s ) in the 
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i-th cell is determined numerically by interval halving; the interval midpoint is taken as the initial 
guess. After the shock location x s is known, the subcell resolution discussed in the previous section 
is applied to the Euler equations in a component-by-component manner. Figure 1 1 illustrates the 
pointwise solution error determined by the modified ENO algorithm on four successively refined 
grids. An odd number of cells has been used for this grid-refinement study in order to facilitate 
the subcell resolution process. The shock movement is thereby contained within one cell. (A more 
general approach would require a search algorithm that finds the cell that harbors a shock before 
applying the subcell resolution technique.) Figure 1 1 shows that the solution converges at the same 
rate downstream of the shock as it does upstream, as shown by the L <*, errors in Table VI. 

DISCUSSION 

The foregoing results raise several issues in regard to the development and application of high- 
order-accurate shock-capturing methods. Although the results presented here were obtained with 
a high-order ENO method, they are consistent with those obtained with linear high-order schemes. 
These issues indicate the need for further investigation into the relative merits of high-order-accurate 
shock-capturing schemes. 

The observation that is, perhaps, the most discouraging is the apparent complexity in achieving 
high-order accuracy in the presence of moving shocks. There are at least two significant drawbacks 
to the inclusion of subcell resolution in a numerical algorithm. The first is the added cost. Although 
the additional expense for the subcell resolution was relatively minor for the sound-shock interaction 
problem, this modification does not suffice for the more general case. For example, the existence 
and location of a single shock was assumed to be known a priori. This assumption was valid 
because of the assuredness that the shock would not move outside the initial cell location, within 
the required sound-wave amplitude range and mesh-refinement parameterization. However, in 
the more general case, every cell must be tested for the harboring of a discontinuity before the 
subcell resolution can be applied. (Such a test is proposed in Ref. 21.) The second and more 
significant drawback to this approach is that its extension to multiple dimensions is not straight 
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forward. A shock is a curve in two-dimensional space, and a surface in three dimensions. The 
multiple parameterization required to extend Eq. 9 to two- and three-dimensional problems would 
be cumbersome, to say the least, and would still not guarantee a unique solution for the shock 
location and shape. 

These observations raise considerable concern in regard to the use of high-order-accurate 
methods in the study of unsteady flows with shocks. If high-order methods only yield first-order 
results, why use them? Before this question can be fully answered, it must first be determined 
whether the first-order error from a high-order method is significantly smaller than that of a lower 
order method. Figure 12 depicts the error of the sound-shock interaction problem with a first-order 
upwind method (no subcell resolution). Upon comparison with the Figure 4, it is clear that the 
solution downstream of the shock is more accurate when using the higher-order method. This 
result is by no means conclusive, as there are many other considerations. For instance, what about 
second- or third-order methods? What about the cost of using a given method with respect to its 
accuracy? These and other concerns are topics for future research. 

The ENO methods applied in the present work are designed to capture shocks narrowly without 
oscillations. However, this property alone was found insufficient to produce globally high-order 
results in the solution of an unsteady problem, which suggests that the nonoscillatory properties 
that guarantee the convergence (in the sense of weak solutions) of a nonlinear scheme do not also 
guarantee the rate of that convergence. If this is the case, why not use a linear method? How 
important is the control of a solution’s total variation from the viewpoint of requiring an answer to 
a given problem? Is it sufficient to simply keep a given computation stable in order to extract the 
desired information from a solution? Further research is necessary in regard to the relative merits 
of linear and nonlinear schemes as they are applied to unsteady flows with shocks. Such research 
should include experiments that compare the solution errors of linear and nonlinear schemes with 
respect to their shock resolution, design accuracy, ahd computational efficiency. 
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Table I. Steady-State Entropy Errors 


N c 

II e 111 

ENO-4-3 
r c 1 1 £ | |oo 

r c 

32 

2.736 E-07 


4.154 E-06 


64 

1.800 E-08 

3.93 

4.370 E-07 

3.25 

128 

1.085 E-09 

4.05 

3.671 E-08 

3.57 

256 

6.372 E-ll 

4.09 

2.697 E-09 

3.77 


Table IV. Solution Errors for IBVP 3 


N c 

ii e Ik 

ENO-4-3 

r c 1 1 e | |oo 

r c 

65 

2.924 E-05 


9.217 E-05 


129 

1.555 E-05 

0.92 

4.825 E-05 

0.94 

257 

7.821 E-06 

1.00 

2.440 E-05 

0.99 

513 

3.910 E-06 

1.00 

1.224 E-05 

1.00 


Table II. Leo Pressure Errors for IBVP 2 Table V. Solution Errors for IBVP 3 


No 

x < 0.45 

ENO-4-3 
r c x > 0.55 

r c 

No 

ENO-4-3-SR 

II e 111 l-c II e lloo 

r c 

64 

8.358 E-05 


1.677 E-03 


65 

2.599 E-06 


1.153 E-05 


128 

6.540 E-06 

3.68 

1.392 E-04 

3.59 

129 

1.552 E-07 

4.11 

8.644 E-07 

3.78 

256 

4.758 E-07 

3.78 

3.087 E-05 

2.17 

257 

9.772 E-09 

4.01 

5.703 E-08 

3.94 

512 

4.511 E-08 

3.40 

1.689 E-05 

0.87 

513 

6.132 E-10 

4.01 

3.628 E-08 

3.99 


Table III. Solution Errors for IBVP 3 Table VI. L^ Pressure Errors for IBVP 2 




ENO-4-3 




ENO-4-3-SR 



N c 

II e Hi 

r c 

II e lloo 

r c 

No 

x < 0.45 

r c 

a; > 0.55 

r c 


64 

2.464 E-06 


1.460 E-05 


65 

7.407 E-05 


1.853 E-03 



128 

1.665 E-07 

3.89 

9.783 E-07 

3.90 

129 

5.665 E-06 

3.71 

1.831 E-04 

3.34 


256 

1.038 E-08 

4.00 

6.108 E-08 

4.00 

257 

3.623 E-07 

3.97 

1.288 E-05 

3.83 


512 

6.383 E-10 

4.02 

3.759 E-09 

4.02 

513 

2.667 E-08 

3.76 

8.266 E-07 

3.96 
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FIGURES 
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Figure 1: a) Nozzle area, b) Mach number, ENO-4-3. 
c) S = P/p 7 , ENO-4-3. 


Figure 3: Solution of the sound-shock interaction problem at 
t = 30 T x , ENO-4-3. 
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Figure 4: Pointwise error for the sound-shock interaction 
problem, ENO-4-3. 
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Figure 2: Quasi- ID Nozzle: Entropy error, ENO-4-3. 


Figure 5: Piecewise polynomial reconstruction. Fluxes are 
determined by point values at cell faces. 
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Figure 10: Pointwise error for the IBVP 3 with the disconti- 
nuity within a cell, using ENO-4-3-SR. 
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Figure 11: Pointwise error for the sound-shock interaction 
problem, using ENO-4-3-SR. 
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Figure 12: Pointwise error for the sound-shock interaction 
problem, using a first-order upwind method. 
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ABSTRACT 

The desire to obtain acoustic information from the numerical so- 
lution of a nonlinear system of equations is a demanding proposition 
for a computational algorithm. High-order accuracy is required for 
the propagation of high-frequency, low-amplitude waves. The ac- 
curacy of an algorithm can be compromised by low-order errors that 
naturally occur in the solution of a particular problem. Such errors 
arise from two sources: the presence of discontinuities in the flow 
field or because the geometry on which the problem is defined is not 
eveiy where smooth to the order of the scheme. The performance of 
high-order accuiate essentially non-oscillatory (ENO) schemes on 
piecewise smooth solutions is well documented. Herein, the perfor- 
mance of these methods on smooth solutions defined on piecewise 
smooth geometries is investigated. The propagation of sound in a 
quasi-one-dimensional nozzle is considered as a test case. Some of 
the issues involved in the extension to two spatial dimensions are 
discussed. 

NUMERICAL METHOD 

For the sake of brevity, the necessary details of the ENO schemes 
to be used in this work are presented within the context of a one- 
dimensional scalar equation, 

| u +£/w=° a) 

A control- volume formulation is obtained by integrating Eq. 1 on 
an interval i/ 2 , ^i+ 1 / 2 ] with center X{ and "volume" A X{ . 
The one-dimensional scalar conservation law can then be written 



-1 

Ax* 


/( «(^+l/2> t) ) ~ /( ) 


( 2 ) 


where 

2 r * i + 1/2 

M t ) = -A~r u(x,t)dx (3) 

is the cell average of u on the i-th interval at time t . Temporal 
integration of Eq. 2 is accomplished by one of the Runge-Kutta 
methods of Shu and Osher (1988). The right-hand side of Eq. 2 is 
approximated in a manner similar to that introduced by Harten et 
al. (1987). A brief description follows. 

To approximate the right-hand side of Eq. 2 to high-order ac- 
curacy, the spatial operator must include a high-order pointwise 
approximation to u(x, t) . However, at a given time t , only the 
cell averages in Eq. 3 are available. Therefore, a pointwise “re- 
construction” of the solution from its cell averages is required. To 
this end, let R be an operator which reconstructs the cell averages 
and yields a piecewise polynomial R(x ; u(t)) of degree r — 1 
which approximates to high order, wherever u(x,t) is 

sufficiently smooth. This operator R acts in a piecewise manner in 
that the solution is locally reconstructed within each cell. Let Vi 
denote the polynomial of degree r— 1 which approximates u(x , t) 
in the i-th cell, at time t , i.e. 

Vi(x) = f?(x; , ^2 — 1/2 < x < *£«+ 1/2 (4) 

= u(x,t) -f 0(h r ) 

where h = max* {A#,-}. The specific method used in this work 
is the "reconstruction by primitive" proposedby Harten et al. (1987) 
and is not detailed here. 

This piecewise reconstruction can cause jumps in the approxi- 
mate solution at the cell interfaces that are O (h r ) in smooth regions 
and 0(1) near discontinuities. The fluxes in Eq. 2 are approximated 
by solving the local Riemann problems at the cell interfaces. Thus, 



u t + au x - 0 , a> 0 


r = 4: 
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Figure 1: Preferred reconstruction stencils that result in fluxes 
that are a) one cell upwind and b) one-half cell upwind. 


the right-hand side of Eq. 2 is replaced by its high-order approxi- 
mation, which yields 



-1 
A xi 


fi+l/2(t) ~ ft- 1 / 2(0 


(5a) 


where 


fi+ 1 / 2 ( 7 ) = / Rm ( Vi{x i+ 1 / 2 ) , ^i+i(®i+i/ 2 ) ) (5b) 

and t/#) denotes the flux that is associated with the 

solution of the Riemann problem whose initial states are ur and 
ur . Upon temporal integration with an appropriately high-order 
Runge-Kutta method, the scheme in Eq. 5 is locally r-th-order 
accurate in the L\ sense (Harten et al., 1987). The extensions of 
these schemes to hyperbolic systems and multiple dimensions that 
are discussed in this work were developed by Harten et al. (1987) 
and Casper and Atkins (1993). 

The most unique aspect of the reconstruction operator R is its 
use of adaptive stenciling. That is, the interpolation set used for 
the approximation of u{x^t) within a given cell is allowed to 
shift in an attempt to use the smoothest possible information. In 
this way, ENO schemes can approximate the smooth regions of 
a piecewise continuous function to high-order accuracy without 
the oscillatory behavior that is associated with interpolation across 
steep gradients. Furthermore, adaptive stenciling enables high- 
resolution shock-capturing. Previous research has shown that the 
accuracy of these schemes can degenerate when the stencils are 
allowed to freely adapt (Rogerson and Meiberg, 1990). Further 
research indicates that this accuracy problem can be remedied by 
biasing the stencils toward those that are linearly stable (Shu, 1990; 
Atkins, 1991). For present purposes, the desired reconstruction 
stencils are centered if r is odd and one cell upwind if r is even. In 
this manner, the resulting schemes have an upwind biased flux, as 
shown for the cases r = 4 and r = 5 in Figure 1 . 

It is important to note that the reconstruction-by-primitive oper- 
ator R is not dependent on uniformity of the computational mesh. 


However, a rectangular mesh (or a smooth transformation to a rect- 
angular mesh) is required for its extension to multiple dimensions 
(Casper and Atkins, 1993). Therefore, the ENO methods described 
above have a potentially important role to play in solving problems 
on cartesian meshes. In this work, the application of these methods 
to piecewise smooth geometries is considered. 


ACOUSTIC WAVE IN A 
QUASI-ONE-DIMENSIONAL NOZZLE 

The high-order ENO methods discussed above are now applied 
to the solution of an acoustic wave in a quasi-one-dimensional 
converging-diverging nozzle. The governing equations are the 
quasi-one-dimensional Euler equations: 

^MJ) + ^AF) = H (6a) 

where 
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The variables p,u, P, E, and A are the density, velocity, pressure, 
total specific energy, and nozzle area, respectively. The equation of 
state is 

P = ( 7 - l)p(E- lu 2 ) 

where 7 is the ratio of specific heats which is assumed to have 
a constant value of 1.4. The flow variables are normalized with 
respect to upstream stagnation conditions. The length scale is A# 
with —200 < x < 80. The area distribution is given by: 


M x ) 


Ao , -200 < X < -100 

Bo — Bi cos(i? 2 ‘E) , —100 < x < 19 

DqADxx, 19 <z<80 


where Ao = 134 and the parameters Bk and Dk are determined 
by requiring C 1 continuity at x = —100 and x — 19. The 
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Figure 2: a) Nozzle area A(x), and b) area derivative dA/dx. 
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Figure 3: Initial steady-state solution. 


prescribed area distribution A{x) and its derivative are illustrated 
in Figure 2. This problem is taken from the test cases in (Hardin 
and Ristorcelli, 1995). 

A steady-state solution is obtained by implementing a fourth- 
order (r = 4) ENO scheme until residuals are driven to machine 
zero. The stencils are everywhere biased toward the upwind stencil 
pictured in Figure la, including at the geometric discontinuities. 
Figure 3 illustrates the steady state solution on a mesh of 280 
uniform cells. It should be noted that this numerically converged 
initial condition cannot be obtained with a freely adaptive stencil 
algorithm. 

After the steady state is achieved, an acoustic disturbance is 
introduced at the inlet, x = —200: 

P(t ) = Pi [1 + e sin(w(x/(M + 1) - t ) )] 

P(*)] 1/7 
Pi . 

2 

u(i) = Ui + — -[c(i)-Ci] 

Here the subscript i denotes the steady inlet state, u> is the circular 
frequency, € is the amplitude, and c = y/jP/p is the local sound 
speed. The calculation is performed with w = 0.17T, e = 10" 6 , 
and a Courant number of 0.5 . Time-accurate, nonreflecting numeri- 
cal boundary conditions (Atkins and Casper, 1994) are employed at 
inflow and outflow for both the initial steady solution and the time- 
dependent solution. The inflow is perturbed for 0 < t/T < 30, 
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Figure 5: Entropy error of the acoustically perturbed nozzle, 
using uniform grid spacing. 


where T = 27r/u> is one period of the incoming acoustic wave. 

Figure 4 depicts the pressure perturbation <5P = P(a?,f) — 
P(x, 0) as a function of the nozzle length, at 10 equally spaced 
time intervals, during one period of the incoming acoustic wave, on 
a mesh of 280 cells. Figure 5 shows a grid refinement error study 
for the acoustically perturbed solution on four uniform meshes. The 
error in this figure is determined by deviation from isentropy and 
the plotted variable is <55, where S = P/p 7 . The error is only 
second order near x = —100 and x = 19, as expected, but is 
fourth order away from these points. The lower-order error arising 
from the geometry has remained local. 

The above problem is now reworked, this time using mesh spac- 
ing that is discontinuous at the derivative discontinuities in the 
nozzle. Specifically, for a mesh of N cells, the grid spacing is 
determined by 

f 100/% , -200 < xi < -100 

= l 119 /N 2 , -100 < Xi < 19 (7) 

( 61/JVs , 19 < xi < 80 

where Ni = N/ 2, N 2 = 57V/14, and Ns = N/7. This mesh 
spacing is compared to uniform spacing for N = 280 in Figure 6. 
The staggered mesh is more highly resolved in the upstream region, 



Figure 4: Envelop of pressure perturbation. 


Figure 6: Mesh spacing with 280 cells. 
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Figure 7: Entropy error of the acoustically perturbed nozzle, 
using staggered grid spacing. 

and less resolved downstream than in the uniform case. This stands 
to reason for this problem, as the acoustic wavelength becomes 
longer in the stream wise direction, as the mean flow expands. Fig- 
ure 7 shows the entropy error for the acoustically perturbed solution 
on four meshes determined by Eq. 7. The solution error in this fig- 
ure is qualitatively similar to that on the uniform meshes in Figure 
5. Again, the lower-order error at the geometric discontinuities 
remains local. 

This problem is again worked on the staggered mesh in Eq. 7, 
this time using stencil adaptation near the geometric discontinuities . 
That is, the operator R is not permitted to interpolate across x — 
— 100 or x = 19, and therefore the stencils are forced to be 
one-sided near these points. The entropy error for the acoustically 
perturbed solution in this case is shown in Figure 8. This plot is 
virtually indistinguishable from Figure 7, even though the stencils 
are clearly different, as shown in Figure 9. The pictured stencils are 
those used in the cubic polynomial interpolation of the characteristic 
variable associated with the eigenvalue A = « + c. The variable 
plotted in Figure 9 is the stencil offset j(i) — js(i ), where j(i) 
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Figure 8: Entropy error of the acoustically perturbed nozzle, 
using staggered grid spacing and grid-adaptive stencils. 
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Figure 9: Stencils for the characteristic variable associated with 
A = u + c. 


is the left-most stencil index and js(i) is the left-most index of 
the linearly stable stencil. For the present problem, the value of 
u + c is always positive. Therefore, js (i) = i — 2 as in Figure la. 
A negative value of this offset represents a further upwind-shifted 
stencil, and a positive value represents a downwind-shifted stencil. 

Given the insignificant difference in Figures 7 and 8, one might 
question the need for the geometric stencil adaptation, being that 
it adds to coding complexity while apparently gaining nothing in 
accuracy. This point is well taken for one-dimensional problems. 
However, the ability to geometrically adapt stencils without losing 
accuracy becomes important for multi-dimensional problems, to be 
discussed in the following section. Another question arises as to 
the stability of the scheme when some of the stencils are shifted 
downwind. The stability of a fourth-order upwind scheme with one- 
sided stencils near a boundary has been demonstrated by Atkins and 
Shu (1995). 

TWO-DIMENSIONAL IMPLEMENTATION 

A solution is desired to the scalar hyperbolic equation 

u t + f(u) x + g{u) y = 0 (8) 

Let the set { Cij }, where 

Cij = [$*— i/2i #1+1/2] * [2/7 — 1/2) 2//+1/2] 

denote a rectangular partition of the x-y plane, with (xi, yj) denot- 
ing the centroid of each rectangle. With a semi-discrete formulation 
in mind, a weak solution of Eq. 8 must satisfy 

A X " 

— Uij(t) = — f S-j-1/2 , j CO — fi-l/2 ,j{t) 

+ 9i t j+l/2(t) —9i J-I/ 2 W 
where = A Xi A yj is the area of Cij , and 

l rxi+ 1/2 ryj+ 1/2 

Uij(t) = — / / u{x,y,t)dxdy (10) 

a ij J^i-i/2 ^yj- 1/2 




Figure 10: Discontinuous mesh spacing which requires no grid- Figure 11: Discontinuous mesh spacing which requires grid- 
based stencil adaptation. based stencil adaptation. 


is the cell average of u in Cij at time t. The fluxes are 
ryj+i/2 

/<+i/ 2 .i(*) = / f(u(xi+i/ 2 ,y,t))dy (11a) 

J Vi- 1/2 
[ X i+ 1/2 

9ij+ 1/2 W = / g(u(x y y j+1/2j t))dx (11b) 

''Vi- 1/2 

As in the one-dimensional case, it is necessary to reconstruct 
the solution pointwise from the cell averages given by Eq. 10. 
The two-dimensional extension of the reconstruction-by-primitive 
operator is denoted by R 2 and is a product of two one-dimensional 
operators. On each cell Cij , the solution u(x , y , t ) is approximated 
by a polynomial Vij , 

Vij{x,y) = R?(x,y,u(t)) , x,yeCij 

= R(x-(R(y,U(t))) (12) 

= u(x,y,t) + 0(h r ) 

where h = ma,Xij{Axi, A yj}. 

The fluxes in Eq. 1 1 are computed as follows. The integration in 
Eq. 1 1 is approximated by an appropriate quadrature formula. The 
polynomials { Vij } that are determined by the reconstruction in 
Eq. 1 2 are evaluated at the chosen quadrature points along each cell 
interface. Pointwise fluxes analogous to the one-dimensional flux in 
Eq. 5b are then computed at each quadrature point and then appro- 
priately weighted and summed. The details of this two-dimensional 
method, as well as its extension to systems of conservation laws, 
can be found in (Casper and Atkins, 1993) and are not repeated 
here. 

It is common practice to decompose a computational domain into 
subdomains when solving multi-dimensional problems on complex 


geometries or when efficiency demands that solutions of such prob- 
lems require several levels of local grid resolution. Furthermore, it 
is often inconvenient, if not cost prohibitive, to connect these subdo- 
mains in a manner that is in keeping with a high-order accurate flow 
solver. The nature of such problems has prompted the development 
of more generalized versions of high-order accurate finite- volume 
methods for unstructured geometries (Barth and Fredrickson, 1991; 
Harten and Chakravarthy, 1991; Abgral, 1991). Depending on the 
application, these generalized methods can be quite complex and 
expensive. Of particular concern is an efficient method of stencil 
adaptation, although some recent progress has been made (Suresh 
and Jorgenson, 1995). 

It is, therefore, of present interest to investigate the use of 
high-order methods which can be applied to piecewise rectangular 
meshes. The implementation of the reconstruction operator R 2 
as a product of two one-dimensional operators greatly benefits the 
efficiency of the scheme, relative to the more generalized approach 
required for unstructured meshes. The structured natured of a rect- 
angular mesh allows for more efficient coding. Furthermore, the 
stencil search is more straight forward, as a choice is made in each 
coordinate direction (Casper and Atkins, 1993). 

Three cases of two-dimensional, discontinuous, rectangular 
meshes are considered. The first is one in which the mesh spacing 
A Xi is discontinuous in x, but continuous in y, and A yj is dis- 
continuous in y, but continuous in x, as pictured in Figure 10. In 
this case, there is no need to shift stencils away from x = Xd or 
2 / ^ a s the grid is fully connected. A solution on this mesh 
can be computed to high-order accuracy, disregarding the disconti- 
nuities in A Xi and A yj , analogous to the one-dimensional case in 
the previous section. 

Figure 11 depicts a mesh which differs from that in Figure 10 
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Figure 12: A geometry which requires grid-based stencil adapta- 
tion in addition to generalized reconstruction along the boundary. 

in that A yj is discontinuous in x. In this case, the reconstruction 
operator R 2 cannot use information on both sides of x — Xd 
in a straight forward manner. Therefore, the forcing of one-sided 
stencils near the interface (dashed line) appears to be a more prudent 
application. The flux integrals in Eq. 11a must be distributed over 
the subintervals into which a larger cell face is divided by smaller 
adjoining cells. 

Perhaps the most common type of mesh discontinuity encoun- 
tered when solving problems on rectangular grids is pictured in 
Figure 12. The mesh must abruptly terminate at a curvilinear 
boundary, making cells of non-rectangular shape near the bound- 
ary. In this case, two different reconstruction methods need be 
applied. The rectangular operator R 2 can be applied interior to 
the mesh interface (dashed line). However, near the boundary, the 
more generalized reconstruction (Harten and Chakravarthy, 1991) 
must by implemented, as the mesh is unstructured in this region. 
Again, the fluxes must be appropriately distributed amongst adjoin- 
ing cells. Furthermore, the fluxes along the boundary must take 
into account the boundary’s curvature (Casper and Atkins, 1993). 

CONCLUDING REMARKS 

ENO methods that employ reconstruction-by-primitive operators 
can be useful in solving computational aeroacoustics problems on 
discontinuous meshes. The retention of accuracy has been demon- 
strated by the fourth-order results obtained in the solution of an 
acoustic wave in a quasi-one-dimensional nozzle. Suggestions have 
been made for the multi -dimensional application of these meth- 
ods to piecewise smooth rectangular meshes. Multi-dimensional 
test cases must be investigated, however, in order to support these 
claims. 


ACKNOWLEDGEMENT 

The author acknowledges the support of NASA Langley Re- 
search Center in Hampton, Virginia, under Grant NAG1-1653. 

REFERENCES 

Abgrall, R., 1991, “Design of an Essentially Non-Os dilatory Re- 
construction Procedure on Finite-Element Type Meshes,” NASA 
Contractor Report 189574, ICASE Report No. 91-84. 

Atkins, H., 1991, “High-Order ENO Methods for the Unsteady 
Navier-Stokes Equations,” AIAA Paper 91-1557. 

Atkins, H. and Casper, J., 1994, “Non-Reflective Boundary Condi- 
tions for High-Order Methods,” AIAA Journal , Volume 32, No. 3, 
pp. 512-518. 

Atkins, H. and Shu, C. W., 1995 “GKS and Eigenvalue Stability 
Analysis of High-Order Upwind Schemes,” in preparation. 

Barth, T. J. and Frederickson, P. O., 1990, “High-Order Solution 
of the Euler Equations on Unstructured Grids Using Quadratic 
Reconstruction,” AIAA Paper 90-0013. 

Casper, J. and Atkins, H. “A Finite- Volume High-Order ENO 
Scheme for Two-Dimensional Hyperbolic Systems,” 1993, Jour- 
nal of Computational Physics , Vol. 106, pp. 62-76. 

Hardin, J. C. and Ristorcelli, J. R., eds., 1994, Proceedings, 
ICASE/LaRC Workshop on Benchmark Problems in Computational 
Aeroacoustics , Hampton, Virginia. 

Harten, A., Engquist, B., Osher, S., and Chakravarthy, S., 
1987, “Uniformly High Order Accurate Essentially Non-oscillatory 
Schemes III,” Journal of Computational Physics , Vol. 71, No. 2, 
pp. 231-323. 

Harten, A. and Chakravarthy, S., 1991, “Multi-Dimensional ENO 
Schemes for General Geometries,” NASA Contractor Report 
187637, ICASE Report No. 91-76. 

Rogerson, A. and Meiberg, E., 1990, “A Numerical Study of the 
Convergence Properties of ENO Schemes,” Journal of Scientific 
Computing, Vol. 5, No. 2, pp. 151-167. 

Shu, C., “Numerical Experiments on the Accuracy of ENO and 
Modified ENO Schemes,” 1990, Journal of Scientific Computing , 
Vol. 5, No. 2, pp. 127-150. 

Shu, C. and Osher, S., “Efficient Implementation of Essentially 
Non-Oscillatory Shock-Capturing Schemes,” 1988, Journal of 
Computational Physics , Vol. 77, No. 2, pp. 439-471. 

Suresh, A. and Jorgenson, P. C. E. 1995, “Essentially Nonoscilla- 
tory (ENO) REcons tractions via Extrapolation,” AIAA Paper 95- 
0467. 


6 



COMPUTATIONAL CONSIDERATIONS FOR THE SIMULATION 
OF DISCONTINUOUS FLOWS 


MARK H. CARPENTER 
NASA Langley Research Center 
Hampton , Virginia 

AND 

JAY H. CASPER 

Old Dominion University 

Norfolk f Virginia 


Abstract. 

The numerical study of aeroacoustic problems places stringent demands 
on the choice of a computational algorithm, because it requires the ability 
to propagate disturbances of small amplitude and short wavelength. The 
demands are particularly high when shock waves are involved, because the 
chosen algorithm must also resolve discontinuities in the solution. In a pre- 
vious work [1] the capabilities and deficiencies of shock-capturing methods 
for aeroacoustic problems were demonstrated using a high-order essentially 
nonoscillatory (ENO) numerical method. It was shown that first-order re- 
sults are obtained when simulating time-dependent flows with discontinu- 
ities. The present study reaffirms this conclusion by comparing the ENO 
results with those obtained using a conventional linear scheme. 

A sixth-order-accurate compact implicit finite difference scheme is used 
to investigate various discontinuous flows. The design order of accuracy 
is achieved in the smooth regions of a steady-state, quasi-one-dimensional 
Euler test case, as well as in the time-dependent Burgers’ equation. How- 
ever, in the unsteady Euler sound-shock interaction, first-order results are 
obtained downstream of the shock. A comparison is made between the 
linear and nonlinear results, noting the advantages of each method. A dis- 
continuous linear model problem is then used to identify the cause of the 
first-order results. Here, the nature of the solution error is quantified as 
being predominantly a numerical phase shift, and a post-processing proce- 
dure is demonstrated which increases the solution accuracy downstream of 
the discontinuity to second-order. 
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1. Introduction 

This work is motivated by the desire to develop numerical methods that 
will be useful in the study of aeroacoustic phenomena that occur in flows 
with shocks. For example, shocks in jet flows, on wings, and in supersonic 
combustion inlets contribute significantly to sound generation. Problems 
such as these represent some of the more challenging aspects of ongoing re- 
search in the developing area of computational aeroacoustics (CAA) . For a 
computational algorithm, obtaining acoustic information from a numerical 
solution that involves shock waves is a demanding proposition. In general, 
high-order accuracy is required for the propagation of high-frequency low- 
amplitude waves. In addition, the shock must be adequately captured. 

In reference [1], a discussion was opened on the relative merits of the 
numerical methods which can simulate sound sources that are generated in 
flows with shocks. A fourth-order essentially nonoscillatory (ENO) numeri- 
cal scheme was used to demonstrate the nature of solution error for steady 
and unsteady cases. It was shown that the quasi-one-dimensional Euler 
equations, in the presence of a shock, admit a steady-state solution which 
is of design order-of-accuracy. The time-dependent Euler equations revert 
to first-order accuracy downstream of the shock. The magnitude of the 
error, however, is much larger when using a uniformly first-order scheme. 

Given the apparent inability of nonlinear schemes to recover design 
accuracy downstream of an unsteady shock, one must reassess the potential 
benefits of nonlinear schemes relative to linear schemes. In this work we 
begin to open the discussion on the relative merits of linear and nonlinear 
schemes for discontinuous flows. In addition, we identify a larger class of 
problems for which design accuracy can be obtained. Finally, we quantify 
the nature of the solution error downstream of a discontinuity when design 
accuracy is not obtained. 

In the following sections, we introduce a fourth-order-accurate essen- 
tially nonoscillatory (ENO) method and a sixth-order compact method, 
and then use them to simulate a steady-state nozzle flow with a shock. In 
both cases design accuracy is obtained away from the discontinuity. It is 
then shown, using Burgers’ equation, that design accuracy can be achieved 
in the time-dependent case, even for nonlinear equations. We then show that 
accuracy suffers if either numerical scheme is used on the time-dependent 
sound-shock interaction problem. The disappointing results in regard to 
the accuracy of this solution are explained through the study of a simpler 
linear model problem. The solution error is quantified as being phase error, 
to leading truncation order. We demonstrate in the linear case that post- 
processing the numerical solution, with the explicit knowledge of the phase 
error, increases the solution error to second order. Some of the difficult is- 
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sues and ramifications for current methods that are raised by these results 
are discussed in the final section. 


2. High-Order Accurate Shock-Capturing 

Many methods are available in the literature that attempt to balance the 
properties ,of high-order accuracy and shock capturing. They can be clas- 
sified into two categories : linear and nonlinear. We will not attempt to 
describe the details of either class of numerical methods, and refer those 
who are interested to a more complete discussion in reference [1]. Rather, we 
present only a brief theoretical review of linear and the nonlinear strategies. 

Within the linear class of numerical shock-capturing schemes, the inter- 
polation set for the approximation of the solution or its derivatives is fixed 
as a function of grid location. Linear methods admit strong oscillations 
in regions in which physical gradients are inadequately resolved. Central- 
differencing operators and spectral methods are particularly prone to these 
numerical oscillations. For problems with discontinuities, limiters or filters 
are usually required to keep oscillations from growing without bound. 

In the nonlinear class of schemes, the strategy with respect to discon- 
tinuities is to employ some sort of adaptive interpolation. The goal is to 
achieve formal high-order accuracy in smooth regions and high shock res- 
olution without oscillations. The class of ENO schemes [2] - [3] has been 
designed to have such properties. As originally presented, the local polyno- 
mial approximation operator adapts its interpolation set to the smoothest 
available part of the solution. 

Although discontinuous solutions generated by a linear strategy are 
usually not as pictorially pleasing as solutions in which shock profiles are 
monotone, linear schemes are more computationally efficient than nonlin- 
ear schemes. The efficiency of a numerical algorithm is extremely important 
for aeroacoustic simulations because such problems are time dependent and 
require a fine computational mesh for the resolution of high-frequency dis- 
turbances. Because nonlinear methods are designed to avoid the production 
of spurious oscillations, the stability of a calculation of a flow with shocks 
is more readily obtained. However, their adaptive interpolation operator 
significantly hampers their efficiency relative to linear schemes. 


3. Steady Shock in a Nozzle 

A steady-state flow with a shock in a quasi-one-dimensional converging- 
diverging nozzle is numerically investigated. The governing equations are 
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the quasi-one-dimensional Euler equations: 

(la) 

where 

' P ' 

U — pu , F= 

.PE. 

The variables p, n, P, P, and A are the density, velocity, pressure, total 
specific energy, and nozzle area, respectively. The equation of state is 

P={ 1 -l)p{E- l -u 2 ) 

where 7 is the ratio of specific heats, which is assumed to have a constant 
value of 1.4. 

The spatial domain of the nozzle is 0 < x < 1 . The nozzle shape is 
determined exactly through the requirement of a linear distribution of Mach 
number from M = 0.8 at the inlet to M = 1.8 at the exit, assuming the 
flow is isentropic and fully expanded. A schematic of the area distribution 
can be found elsewhere [1]. 

Given the prescribed area distribution, the Mach 0.8 inflow state is 
retained at # = 0, and the outflow condition at x = 1 is determined such that 
a shock forms at x s = 0.5, which corresponds to a preshock Mach number 
of M = 1.3 . Steady-state solutions are obtained by implementing a fourth- 
order accurate finite- volume ENO scheme and a sixth-order finite-difference 
scheme until residuals are driven to machine zero. In the ENO algorithm, 
the spatial accuracy is achieved by solving the equations in control-volume 
form as presented in Ref. [2]. The equations are integrated in time via a 
third-order-accurate Runge-Kutta scheme [3]. This numerical method will 
be referred to as “ENO-4-3.” As has been established in previous research, 
[4] - [7] the adaptive stencils employed in the spatial operator are biased in 
smooth regions toward those that are linearly stable. 

The finite-difference algorithm uses a compact implicit spatial opera- 
tor, shown elsewhere to be time-stable and spatially a sixth-order accurate 
discretization [8]. A small amount of background tenth-order numerical 
dissipation was added to stabilize the calculation. No explicit filtering is 
needed in these calculations. The time advancement scheme is a five stage 
fourth-order low-storage RK scheme [9]. This scheme which is formally 
sixth-order in space and fourth-order in time will be referred to as the 
“LIN-6-4” scheme. 

One of the simpler methods of determining the error of this solution 
relies on the fact that the value of the entropy-like quantity S = P/p 1 is 
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Figure 1. Quasi-lD Nozzle: Entropy error, ENO-4-3. 


piecewise constant. The exact magnitude of the jump is known analytically. 
The pointwise entropy error for this solution on four successively refined 
meshes is illustrated in Fig. 1. It is evident from Fig. 1 that the accuracy is 
fourth order on either side of the shock (error decay by a factor of 16 with 
each grid doubling). 



Figure 2. Quasi- ID Nozzle: Entropy error, LIN-6-4, 

Fig. 2 shows the steady-state Entropy error as determined by the LIN- 
6-4 scheme. Note that an order one error always exists at the location of 
the discontinuity. This error which is not present in the ENO-4-3 results, 
is inherent in the linear scheme as a consequence of the Gibbs’ oscillations. 
Away from the discontinuity the scheme converges at the design accuracy 
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of sixth-order. The ENO-4-3 scheme is considerably more accurate than the 
LIN-6-4 scheme on coarse grids, despite its lower formal accuracy. Similar 
results were obtained for test cases in which the shock location did not 
coincide with a grid point. 

The time independence of the solution makes this a convenient example 
for the demonstration of high-order accuracy in the presence of a shock. 
Both the linear and nonlinear schemes were able to obtain higher-order 
accuracy away from the discontinuity, and the nonlinear schemes’ behavior 
in the neighborhood of the shock was exceptionally good. We assert that 
these results are representative of linear and nonlinear schemes, and are 
general for the steady-state Euler equations in one spatial dimension. (The 
generalization to multiple spatial dimensions is not straightforward and is 
the topic of a forth-coming paper). 

4. Time Dependent Interactions 

A moving shock presents a greater challenge in regard to high-order-accurate 
shock capturing. As will be seen, many of the desirable attributes of the 
nonlinear schemes are lost in the time-dependent case. We begin our discus- 
sion with a simple test case in which the equation is nonlinear and supports 
a time-dependent discontinuity. 

4.1. BURGERS’ EQUATION 

Studying Burgers’ equation helps to clarify an important issue about solu- 
tion accuracy in the presence of discontinuities. It is common to use Burgers’ 
equation as a model problem to mimic the nonlinear nature of the Euler 
equations. The quadratic nonlinearity of the spatial flux in Burgers’ equa- 
tion is similar to the convective terms in the Euler equations. An assertion 
frequently implied in the literature is that solution error will be qualita- 
tively similar for Burgers’ equation and the Euler equations. We show by 
counter-example that the nature of the solution error for discontinuous 
flows is fundamentally different. Thus, Burgers’ equations is of only limited 
value in studying time-dependent discontinuous flows. 

In the absence of viscosity, Burgers’ equation can be expressed as ^ + 
= o. We employ initial conditions U(x,0) = 1; 0 < x < [/(a;,0) = 

— 1; - < x < 1 and boundary conditions 17(0, t) = 1 + esin(uit ); t > 0, 
[7(1, t) = —1 + esin{u 2 t)\ t > 0. Here e = 0.01, = 127T, and o >2 

is related to by a phase shift. A simulation using the LIN-6-4 scheme 
is run to a physical time T = 10. An “exact” solution is obtained using 
a two-domain Chebyshev discontinuity fitting code. Sufficient resolution in 
space and time is used to ensure that further refinement does not change 
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the interpolated solution on the uniform grid. 

Large oscillations near the discontinuity can be observed in the numer- 
ical solution obtained with the LIN-6-4 scheme. Fig. 3 shows the pointwise 
error on four successive grids, plotted on a Logarithmic scale. Qualitatively, 
the solution error is similar to that obtained in the steady-state nozzle us- 
ing the LIN-6-4 scheme. An order one error exists at the discontinuity, but 
the solution converges at the design rate of sixth-order on both sides of 
the discontinuity. We note here, and clarify later in this section, that this 
high-order behavior is not achieved in more general hyperbolic equations. 



Figure 3. Burgers’ equation: solution error, LIN-6-4. 


4.2. SOUND-SHOCK INTERACTION 

The effects of shocks on sound waves, (and vice versa), are important to 
the acoustics and performance of aircraft design. Therefore, the ability to 
obtain an accurate solution to such a model initial-boundary-value problem 
(IBVP) is important in the development of shock-capturing methods for 
CAA research. 

The governing equations are the one-dimensional Euler equations: 

^ + A F(£0 = o (2a) 

where the components of U and F(U) are identical to those given in Eq. 
lb. The equation of state is also the same as in the previous example. 

The spatial domain is 0 < x < 1 . The piecewise constant initial con- 
ditions, U L and Ur , are those of a steady shock located at x s = 0.5. The 
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flow is from left to right, and the state Ul is a Mach 2 flow upstream of 
the shock. The flow variables are normalized with respect to this upstream 
flow. At t = 0, an acoustic disturbance is introduced at x = 0: 


p(o,t) 

u( 0, t) 


= Pl ( 1 + e sin cot ) 

rp(o,t)i 1/7 


= PL 


UL 


Pl 

2 


+ ^Tj;[ c ( 0, f ) _ CL ] 


(2b) 


where u is the circular frequency, e is the amplitude, and c = \fyPjp is the 
local sound speed. 

The numerical solution of this problem is obtained through the imple- 
mentation of both the ENO-4-3 and the LIN-6-4 algorithms. Both algo- 
rithms should achieve their design spatial accuracy for a suitably small 
temporal step size. The exact solution is obtained by a two-domain Cheby- 
shev spectral technique [9]. Shock fitting is used to divide the domain into 
two computational regions. A Chebyshev collocation method is used in each 
region for the spatial discretization. A fourth-order Runge-Kutta scheme is 
used to discretize time. Sufficient spatial and temporal resolution are used 
to guarantee machine precision of the solution. 


5 P 



Figure ]±. Solution of the sound-shock interaction problem at t = 30 7\, ENO-4-3. 

Fig. 4 depicts the pressure perturbation SP(t ) = P(x , t) — P(x , 0) at t = 
30 Tx, where T\ = 2k / u is one period of the incoming acoustic wave. The 
acoustic wave amplitude is e = 0.001, and u = 27 rk(uL + cl) is determined 
by requiring a wave number k = 6 with respect to unit length and a mean 
wave speed u^ + cl. The calculation, represented by circles, was performed 
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on a uniform mesh of 256 cells and a Courant number of 0.5, with the ENO- 
4-3 code. The exact solution is represented by a continuous line. In this 
pictorial measure, the numerical algorithm performs well with respect to 
its prediction of the amplified sound wave at higher frequency downstream 
of the shock. The missing circle values near the shock are off the plot and are 
due to the use of the stencil-biasing parameters near a moving discontinuity. 


log I Error I 



Figure 5. Pointwise error for the sound-shock interaction problem, ENO-4-3. 


Even more instructive, however, is the pointwise error made by this 
calculation with respect to the mesh width. Fig. 5 illustrates this error 
on four successively refined meshes. The solution is clearly fourth-order 
accurate upstream of the shock, but only first-order downstream of the 
shock, as shown by the L <*, error data in Table 1. The errors are computed 
on two spatial subdomains: 0 < x < 0.45 and 0.55 < x < 1. In this manner, 
the first-order error that is generated in the vicinity of the shock is avoided. 


TABLE 1. Loo Pressure Errors : ENO-4-3 


N c 

a; <0.45 

r c 

x >0.55 

r c 

64 

8.358 E-05 


1.677 E-03 


128 

6.540 E-06 

3.68 

1.392 E-04 

3.59 

256 

4.758 E-07 

3.78 

3.087 E-05 

2.17 

512 

4.511 E-08 

3.40 

1.689 E-05 

0.87 


Fig. 6 depicts the solution error for the sound-shock interaction problem, 
using the LIN-6-4 scheme. (Here, in addition to the tenth-order damping, 
the solution was filtered with a tenth-order filter every 5 time-steps). The 
qualitative features of the solution and its error are nearly identical to 
those obtained with the ENO-4-3 scheme. The order one error near the 
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discontinuity is a bit more localized in the ENO-4-3 case, but the magnitude 
of solution error downstream of the discontinuity is nearly the same. The 
design accuracy is obtained upstream of the discontinuity, but the solution 
accuracy is first-order downstream of the shock. 



Figure 6. Point wise error for the sound-shock interaction problem, LIN-6-4. 

Comparing the error differences in the Euler and Burgers’ equations can 
lend insight into when a first-order solution will occur at a time-dependent 
discontinuity. We begin by noting how error propagates in each respective 
equation. In Burgers’ equations, all information propagates towards the 
discontinuity. Thus, error which is generated at the discontinuity is localized 
at the discontinuity. In the Euler equations, error that is generated at the 
discontinuity is swept downstream on the down running characteristics. (It 
can be shown that all information traveling towards the discontinuity is 
design-order accurate on both sides of the discontinuity, until nonlinear 
corruption occurs near the shock). Information traveling away from the 
discontinuity is first-order. All flow variables being combinations of the 
two, are no better than first-order accurate. Burgers’ equation is not a 
good model problem for determining the error characteristics of a time- 
dependent numerical scheme. 

The general observation that both linear and nonlinear schemes are 
first-order downstream of the shock raises considerable concern in regard 
to the use of high-order-accurate methods in the study of unsteady flows 
with shocks. If high-order methods only yield first-order results, why use 
them? Before this question can be fully answered, it must first be deter- 
mined whether the first-order error from a high-order method is signifi- 
cantly smaller than that of a lower order method. Fig. 7 depicts the error 
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of the sound-shock interaction problem with a first-order upwind method. 
Upon comparison with the Fig. 5 or Fig. 6, it is clear that the solution 
downstream of the shock is more accurate when using either higher-order 
method. This result is by no means conclusive, as there are many other 
considerations. For instance, what about second- or third-order methods? 
What about the cost of using a given method with respect to its accuracy? 
These and a more general comparison of linear and nonlinear scheme are 
topics for future research. 

log 1 Error I 



Figure 7. Pointwise error for the sound-shock interaction problem, using a first-order 
upwind method. 

The first-order nature of the solution downstream of the discontinuity is 
difficult to study using the full Euler equations. In the next section, we con- 
struct a simple model problem which contains the underlying mechanism 
for the degradation to first-order accuracy. 

5* A Linear Model Problem 

The lower-order accurate results of the previous section can be analyzed 
through the study of a discontinuous linear scalar model problem, wherein 
we isolate the important phenomenon of propagation of information through 
a discontinuity. This trait will be common to almost any aeroacoustic prob- 
lem that involves shock waves. Consider the scalar equation: = 

0, where the piecewise constant wave speed a(x) is a(x) = , x < | and 

a(x) = or , x > - on the interval 0 < x < 1. (For the finite difference 
code, the jump value at a; = | is assigned to be the average of the left 
and right states.) The initial condition is given by w(ar,0) = sin(a;x) on 
0 < cc < 1, and the boundary condition is w(0,i) — sin(u> (—ait) at 
x = 0, t > 0. The constant a ; is determined by requiring a wave number 
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k = 4 with respect to unit length and the downstream wave speed clr. The 
exact solution is given at long times (times after the initial solution hats 
swept out of the domain) by u{x,t) — sinfa; (x — a^t)] on the interval 
0 < x < I an ^ u ( x it) = f£sin[o; (-£(& — |) — a£,£ + §] on the interval 
- < a; < 1. The exact solution is obtained by requiring conservation of 
fluxes at the interface. 

The loss in accuracy in numerical solutions of linear problems with 
discontinuous initial data has been the subject of previous research by other 
authors [11] - [14]. All of these previous studies involved solutions of coupled 
linear systems. It is, therefore, instructive to note that the solution of this 
discontinuous linear problem is analogous to that of a coupled system in the 
following way. The flux conservation condition couples at the discontinuity 
location the equations in the left and right domains. The coupling occurs 
at a discrete point (like a boundary) and the information transfer is from 
one equation to the other. As in the present study, the work of Majda and 
Osher [11] was concerned with the inherent degradation in accuracy in a 
region where information is numerically propagated across a discontinuity. 



Figure 8. Point wise error for discontinuous linear model problem run with LIN-6-4 

Fig. 8 shows the pointwise error for the discontinuous linear model prob- 
lem with a wave-speed ratio of ^ = 2. The simulation is run to time T = 4 
on a series of six grids with the £ IN-6-4 scheme. The solution is filtered ev- 
ery five time steps. As with the Euler equations, the solution accuracy is 
sixth-order upsteam, and first-order downstream of the discontinuity. 

Inspection of the numerical solution downstream of the discontinuity re- 
veals that the dominant component of the error is related to phase. Specifi- 
cally, the amplitude of the waves downstream of the discontinuity is nearly 
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correct, but the waves are displaced by some fraction of the grid spacing. 
To quantify this behavior, the Fourier transform of the solution is formed 
on the interval | < x < 1. The exact transform of the downstream sine 
wave on the interval (made periodic by construction) is the k = 1 mode of 
unit amplitude. Fig. 9 plots the amplitude of the Fourier modes obtained 
from the numerical solution on this interval. The energy in the spurious 
higher modes decays very rapidly, indicating that the amplitude portion of 
the solution error is higher-order. 



log k 


Figure 9, Solution error as a function of Fourier wave number k 

The phase shift of the numerical solution as a percentage of grid spacing 
(as determined by the phase angle of the Fourier transform), asymptotes 
to a nonzero constant, however. For sufficient resolution, the numerical 
solution is shifted a finite percentage of the grid spacing, independent of 
the resolution. If the numerical solution is adjusted to account for this phase 
shift, then the solution error can be reduced. Fig. 10 plots the data shown 
in Fig. 8 where the exact solution downstream of the discontinuity is shifted 
by a constant multiple of the grid spacing. It is apparent that the solution 
accuracy is increased by this procedure. Unfortunately, the convergence 
rate downstream of the discontinuity is still only second-order and not the 
design order of six as would be expected with the LIN-6-4 scheme. 

A more troubling observation about the numerical phase shift, is that 
the shift is amplitude dependent. Table 2. shows the dependence of the 
phase shift on the wave-speed ratio Multiple waves interact at a shock 
in the Euler equations, each having a different wave-speed ratio. To correct 
the multiple wave situation would require component by component shifting 
and would be more difficult to implement. 
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Figure 10. Pointwise solution error after correcting for the numerical phase shift, with 
LIN-6-4 


TABLE 2. Phase shift verses 


u L 

6 

Ur 

Ax 

4/1 

0.2351 

2/1 

0.08412 

1/2 

0.1700 

1/4 

0.9674 


6. Discussion 

In this work, and elsewhere [1], we provide a systematic study of the ef- 
fects of discontinuities on solution accuracy. We focus on the ENO-4-3 and 
LIN-6-4 formulations as being representative of high-order nonlinear and 
linear schemes in use today. We show for higher-order formulations, that 
the solution accuracy is strongly dependent on the mathematical character 
of the governing equations, and only weakly dependent on the numerical 
method. 

Specific conclusions include: 1) The quasi-one-dimensional Euler equa- 
tions admit design accuracy away from a steady discontinuity. The ENO-4-3 
results are significantly better than the LIN-6-4 results for the steady case. 
2) The unsteady Burgers’ equation admits design accuracy away from a 
discontinuity, in contrast to the unsteady Euler equations. Equations that 
focus information at the discontinuity maintain design accuracy away from 
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that discontinuity. 3) The unsteady Euler equations are first-order accurate 
downstream of a discontinuity. The ENO-4-3 and LIN-6-4 schemes achieve 
similar solution accuracy downstream of the shock for the sound-shock in- 
teraction problem. 

We use a linear model problem with a discontinuous wave speed, to 
study the nature of discontinuity related solution error. We recover first 
order solution error downstream of a discontinuity in the linear case, and 
quantify the error as being phase related to leading truncation order. Second- 
order results, are obtained by post-processing the solution with the known 
phase shift. It is doubtful whether this post-processing step is general for 
systems of equations in multiple dimensions. 

Research continues in several directions. It is questionable whether de- 
sign accuracy can be obtained for the discontinuous steady Euler equa- 
tions in multiple spatial dimensions. (The mathematical character of the 
2-D steady and the 1-D unsteady Euler equations is similar, and we know 
the later to be first-order downstream of discontinuities). Further inves- 
tigation continues into the relative merits of high-order-accurate shock- 
capturing schemes, and in particular, the relative merits of linear and non- 
linear schemes. Attention focuses on shock resolution, design accuracy, and 
computational efficiency. Work continues into developing general techniques 
which recover design accuracy in the presence of discontinuities. Linear 
schemes suffer from a first order phase shift which is wave-speed jump de- 
pendent. It is possible, therefore, to correct the time-dependent solution for 
the full Euler equations making use of phase shift information. 

We conclude by noting that the inability to recover design accuracy 
downstream of a discontinuity is not a proof of inherent first-order accuracy; 
but rather only a demonstration that design accuracy has not been found! 
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